%%time
import datetime
import os, glob, sys
import warnings
warnings.filterwarnings('ignore', '.*invalid value encountered in true_divide.*', )
import pickle
import numpy as np
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
from matplotlib.colors import ListedColormap
from matplotlib import cm
from scipy.interpolate import UnivariateSpline
from sklearn.metrics import r2_score, mean_squared_error, explained_variance_score
from scipy.stats import pearsonr, spearmanr
mpl.rcParams['savefig.dpi'] = 300
CPU times: user 1.39 s, sys: 1.96 s, total: 3.35 s Wall time: 1.57 s
%%time
print(datetime.datetime.now())
dirroot = '/mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/'
dirshared = dirroot + 'SHARED-DATAS/'
dircmip6 = dirroot + 'CMIP6-DATAS-LINKS/'
2023-11-28 17:15:04.429097 CPU times: user 19 µs, sys: 21 µs, total: 40 µs Wall time: 41 µs
%%time
print(datetime.datetime.now())
dirout = dirroot + 'dimred-220708-figure-article-v0/'
if not os.path.isdir(dirout) : os.mkdir(dirout)
2023-11-28 17:15:04.434420 CPU times: user 207 µs, sys: 0 ns, total: 207 µs Wall time: 174 µs
%%time
print(datetime.datetime.now())
sys.stdout.echo = open(dirout+'stdout.txt', 'w')
sys.stderr.echo = open(dirout+'stderr.txt', 'w')
2023-11-28 17:15:04.439044 CPU times: user 0 ns, sys: 639 µs, total: 639 µs Wall time: 46.6 ms
mpl.rcParams['font.size'] = 8
cm2in = 1/2.54
%%time
print(datetime.datetime.now())
print('Maps score KR: inputs')
data2plot_dir = dirroot + 'dimred-220602-score-analysis/'
grid = '360x180'
data2plot_file = data2plot_dir+'data2plot-maps-score-KR-'+grid+'.pckl'
# model_list = ['NorESM2-LM', 'MPI-ESM1-2-LR', 'CNRM-ESM2-1', 'CESM2', 'ACCESS-ESM1-5', 'IPSL-CM6A-LR']
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
#grid_list = ['90x45', '180x90']
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2022-12-05 18:01:24.738793 Maps score KR: inputs CPU times: user 349 µs, sys: 0 ns, total: 349 µs Wall time: 218 µs
%%time
print(datetime.datetime.now())
print('Maps score KR: fig')
#-----------------
# Load data2plot
#-----------------
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (18*cm2in, 16*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
from matplotlib.colors import ListedColormap
from matplotlib import cm
cmap = cm.get_cmap('viridis', 256)
newcolors = cmap(np.linspace(0, 1, 10))
cmap = ListedColormap(newcolors[:])
cmap.set_under(color='silver')
cmap.set_over(color='silver')
cmap.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':1, 'cmap':cmap, \
'transform':ccrs.PlateCarree() }
#---------------------
# Plot score maps
#---------------------
irow = 0
for model in model_list:
icol = 0
for simu in simu_list:
zax = ax[irow, icol]
X = data2plot['maps'][model][simu]['X']
Y = data2plot['maps'][model][simu]['Y']
Z = data2plot['maps'][model][simu]['Z']
pcm = zax.pcolormesh(X, Y, Z, **kwmap)
zax.set_title(subnum.pop()+') '+ model+'\n'+simu_name[simu].upper())
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
icol+=1
#
irow+=1
#
# fig.tight_layout()
#---------------------
# Repositionning axes
#---------------------
for irow, axrow in enumerate(ax[1:]):
for icol, zax in enumerate(axrow):
zw1 = zax.get_position()
ny0 = zw1.y0 - (irow+1)*0.1*zw1.height
zax.set_position([zw1.x0, ny0, zw1.width, zw1.height])
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 1].get_position()
zw2 = ax[0, 2].get_position()
nx0 = zw1.x0
ny0 = zw1.y1 + .9*zw1.height
nw = zw2.x1 - nx0
nh = 0.03*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm, cax=cax, orientation='horizontal', \
ticklocation='top', extend='min')
cbar.set_label(r'Score of the Kernel Regression')
#---------------------
# Save figure
#---------------------
fignam = 'maps-score-KR.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# 23s
2022-12-05 18:01:31.671421 Maps score KR: fig File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220602-score-analysis/data2plot-maps-score-KR-360x180.pckl Figure saved: maps-score-KR.png CPU times: user 27.1 s, sys: 2.82 s, total: 29.9 s Wall time: 22.8 s
%%time
print(datetime.datetime.now())
print('Histogram and interquartile gathering all models for each scenario')
zfact = 1e6*3600*24 # kmol/m2/s -> mmol/m2/d
#-----------------
# FIGURE PARAM
#-----------------
nrow = 2
ncol = 2
fsize = (18*cm2in, 16*cm2in) # width height
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# KEYWORDS DICT
#-----------------
kwhist = dict(bins=50, range=[0, 1], cumulative=True, density=True, color='black')
kwlines = dict(color='gray', ls='--')
#---------------------
# Plot histograms
#---------------------
for isimu, simu in enumerate(simu_list):
zwdata = []
for model in model_list:
zwdata.extend(data2plot['maps'][model][simu]['Z'].flatten())
#
zax = ax.flat[isimu]
pcm = zax.hist(zwdata, **kwhist)
quantile50 = np.nanquantile(zwdata, .5)
quantile20 = np.nanquantile(zwdata, .2)
quantile05 = np.nanquantile(zwdata, .05)
zax.axhline(0.05, **kwlines)
zax.axhline(0.2, **kwlines)
zax.axhline(0.5, **kwlines)
zax.axvline(quantile50, **kwlines)
zax.axvline(quantile20, **kwlines)
zax.axvline(quantile05, **kwlines)
zax.set_title(subnum.pop()+") "+simu.upper()+"\n95pct above: %.2f\n80pct above: %.2f\n50pct above: %.2f\nmean: %.2f"
%(quantile05, quantile20, quantile50, np.nanmean(zwdata)), loc='left')
#
#
fig.suptitle('Histogram and interquartile gathering all models for each scenario', y=1.03)
fig.tight_layout()
2022-12-05 18:02:09.339549 Histogram and interquartile gathering all models for each scenario CPU times: user 3.52 s, sys: 624 ms, total: 4.14 s Wall time: 2.57 s
%%time
print(datetime.datetime.now())
print('Maps summarizing the score of KR: inputs')
data2plot_dir = dirroot + 'dimred-220602-score-analysis/'
grid = '360x180'
data2plot_file = data2plot_dir+'data2plot-maps-summarize-score-KR-'+grid+'.pckl'
# model_list = ['NorESM2-LM', 'MPI-ESM1-2-LR', 'CNRM-ESM2-1', 'CESM2', 'ACCESS-ESM1-5', 'IPSL-CM6A-LR']
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
#grid_list = ['90x45', '180x90']
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2023-07-31 15:25:43.488195 Maps summarizing the score of KR: inputs CPU times: user 73 µs, sys: 3 µs, total: 76 µs Wall time: 80.1 µs
%%time
print(datetime.datetime.now())
print('Maps summarizing the score of KR: fig')
#-----------------
# Load data2plot
#-----------------
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
zfact = 1e6*3600*24 # kmol/m2/s -> mmol/m2/d
gridpt_name = 'gridpt 1'
#-----------------
# FIGURE PARAM
#-----------------
nrow = 2
ncol = 2
fsize = (18*cm2in, 16*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
from matplotlib.colors import ListedColormap
from matplotlib import cm
cmap = cm.get_cmap('viridis', 256)
newcolors = cmap(np.linspace(0, 1, 10))
cmap = ListedColormap(newcolors[:])
cmap.set_under(color='silver')
cmap.set_over(color='silver')
cmap.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
#_________________
# Maps
kwmap = {'vmin':0, 'vmax':1, 'cmap':cmap, \
'transform':ccrs.PlateCarree() }
# kwscatter = dict(edgecolor='k', facecolor='lightgray', alpha = .7, s=40, \
# transform=ccrs.PlateCarree())
kwscatter = dict(color='lightgray', alpha = .6, s=40, \
transform=ccrs.PlateCarree())
kwscatter2 = dict(edgecolor='k', facecolor='none', alpha = 1, s=40, \
transform=ccrs.PlateCarree())
# kwscatter = dict(marker='+', s=60, linewidths=3, \
# transform=ccrs.PlateCarree())
kwrec = {'ls':'--', 'c':'r', 'lw':1, 'transform':ccrs.PlateCarree() }
#_________________
# Time series
kwtruth = dict(ls='-' , label='Target' , color='k')
kwkr = dict(ls='--', dashes=[2, 1], label='Reconstructed' , color='orange')
kwstd = dict(label='$\pm$ std' , color='silver')
kwref = dict(ls='--', label='Trend + Seasonal cycle' , color='royalblue')
kwbbox = dict(boxstyle ="round", fc ="0.8")
#---------------------
# Plot score maps
#---------------------
for isimu, simu in enumerate(simu_list):
zax = ax.flat[isimu]
X = data2plot['maps'][simu]['X']
Y = data2plot['maps'][simu]['Y']
Z = data2plot['maps'][simu]['min score']
pcm = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
# zax.set_title(subnum.pop() + ') ' + simu.upper(), loc='left')
# add regions SO, EP, NA
#from shapely.geometry.polygon import Polygon
#xll, yll, xur, yur = 100, -10, -100, 10
#pgon = Polygon( ( (xll, yll),
# (xll, yur),
# (xur, yur),
# (xur, yll),
# (xll, yll) ) )
#zax.add_geometries([pgon], crs=ccrs.PlateCarree(), facecolor = None, edgecolor='black')
if simu == 'picontrol':
# EP
xll, yll, xru, yru = 125, -10, 280, 10
zax.plot([xll, xru, xru, xll, xll],
[yll, yll, yru, yru, yll],
**kwrec)
# SO
xll, yll, xru, yru = 0, -70, 360, -40
zax.plot([0, 360], [-40, -40], **kwrec)
# NA
xll, yll, xru, yru = -80, 10, -1, 60
xxx = [-65, -13, -5, -1, -50, -78, -65]
yyy = [ 10, 10, 35, 70, 70, 35, 10]
zax.plot(xxx, yyy, **kwrec)
#
# fig.tight_layout()
#---------------------
# Repositionning axes
#---------------------
for icol, zax in enumerate(ax[1, :]):
zw1 = zax.get_position()
ny0 = zw1.y0 - .8*zw1.height
zax.set_position([zw1.x0, ny0, zw1.width, zw1.height])
#
for irow, zax in enumerate(ax[:, 1]):
zw1 = zax.get_position()
nx0 = zw1.x0 - .15*zw1.width
zax.set_position([nx0, zw1.y0, zw1.width, zw1.height])
#
#---------------------
# Plot score timeseries
#---------------------
axts = []
for isimu, simu in enumerate(simu_list):
zw1 = ax.flat[isimu].get_position()
nx0 = zw1.x0+.05*zw1.width
ny0 = zw1.y0 - 1.15*zw1.height
nw = .9*zw1.width
nh = .9*zw1.height
zaxts = fig.add_axes([nx0, ny0, nw, nh])
axts.append(zaxts)
istart = -5*12
X = data2plot['timeseries'][simu][gridpt_name]['time'] [istart:]
Yref = zfact*data2plot['timeseries'][simu][gridpt_name]['Ref'] [istart:]
Ykr = zfact*data2plot['timeseries'][simu][gridpt_name]['KR'] [istart:]
Ykrsup = zfact*data2plot['timeseries'][simu][gridpt_name]['KRsup'][istart:]
Ykrinf = zfact*data2plot['timeseries'][simu][gridpt_name]['KRinf'][istart:]
Ytruth = zfact*data2plot['timeseries'][simu][gridpt_name]['Truth'][istart:]
gridpt = data2plot['timeseries'][simu][gridpt_name]['gridpt']
lltruth, = zaxts.plot(X, Ytruth-Yref, **kwtruth)
#llref, = zaxts.plot(X, Yref , **kwref)
llstd = zaxts.fill_between(X, Ykrinf-Yref, Ykrsup-Yref, **kwstd)
llkr, = zaxts.plot(X, Ykr-Yref, **kwkr)
zaxts.set_xlabel('Years')
# compute pearson correlation
zwYref = zfact*data2plot['timeseries'][simu][gridpt_name]['Ref']
zwYkr = zfact*data2plot['timeseries'][simu][gridpt_name]['KR'] - zwYref
zwYtruth = zfact*data2plot['timeseries'][simu][gridpt_name]['Truth'] - zwYref
zwpearsonr = pearsonr(zwYtruth, zwYkr)
# add min score
jj, ii = gridpt[0], gridpt[1]
zwminscore = data2plot['maps'][simu]['min score'][jj, ii]
bbox = dict(boxstyle="round", fc="0.9")
# zaxts.annotate('score: %.2f\ncorr: %.2f\npval:%.2f' %(zwminscore, zwpearsonr[0], zwpearsonr[1]), (0.65, 0.05),xycoords='axes fraction')
zaxts.annotate('score: %.2f' %(zwminscore), (0.65, 0.05),xycoords='axes fraction')
# Add point on map
jj, ii = gridpt[0], gridpt[1]
X = data2plot['maps'][simu]['X']
Y = data2plot['maps'][simu]['Y']
ax.flat[isimu].scatter(X[jj, ii], Y[jj, ii], **kwscatter)
ax.flat[isimu].scatter(X[jj, ii], Y[jj, ii], **kwscatter2)
#
axts = np.array(axts).reshape((2, 2))
#---------------------
# Set time series labels
#---------------------
for irow, zax in enumerate(axts[:, 1]):
zax.yaxis.set_ticks_position('right')
zax.yaxis.set_label_position('right')
zax.set_ylabel('[mmol$\,$C$\,$m$^{-2}\,$d$^{-1}$]\nCO$_2$ flux')
#
for irow, zax in enumerate(axts[:, 0]):
zax.set_ylabel('CO$_2$ flux\n[mmol$\,$C$\,$m$^{-2}\,$d$^{-1}$]')
#
for isimu, simu in enumerate(simu_list):
modelname = data2plot['timeseries'][simu][gridpt_name]['model name']
zax = ax.flat[isimu]
zaxts = axts.flat[isimu]
zax.set_title(subnum.pop() + ') ' + simu_name[simu].upper())
zaxts.set_title(subnum.pop() + ') %s' %(modelname), loc='left')
#
#---------------------
# Set legend
#---------------------
# zw1 = axts[1, 0].get_position()
# zw2 = axts[1, 1].get_position()
# nx0 = zw1.x1
# nwidth = zw2.x0 - nx0
# nheight= 0.7*zw1.height
# ny0 = zw1.y0 - nheight
# axleg = fig.add_axes([nx0, ny0, nwidth, nheight])
# #axleg.legend(handles=[lltruth, llref, llkr], ncol=3, loc='lower center')
# axleg.legend(handles=[lltruth, llkr], ncol=2, loc='lower center')
# axleg.axis('off')
axts[1, 0].legend(handles=[lltruth, llkr], ncol=1, loc='upper right',
markerfirst=False, handlelength=.7, handletextpad=.3,
frameon=False)
#---------------------
# Colorbar
#---------------------
# zw1 = ax[0, 0].get_position()
# zw2 = ax[0, 1].get_position()
# nx0 = zw1.x0 + 0.3*zw1.width
# ny0 = zw1.y1 + .5*zw1.height
# nw = zw2.x1 - 0.3*zw2.width - nx0
# nh = 0.05*nw
# cax = fig.add_axes([nx0, ny0, nw, nh])
# cbar = fig.colorbar(pcm, cax=cax, orientation='horizontal', \
# ticklocation='top', extend='min')
# cbar.set_label(r'Inter-model minimum score of the Kernel Regression')
zw1 = axts[1, 0].get_position()
zw2 = axts[1, 1].get_position()
nx0 = zw1.x0 + 0.3*zw1.width
ny0 = zw1.y0 - .65*zw1.height
nw = zw2.x1 - 0.3*zw2.width - nx0
nh = 0.05*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm, cax=cax, orientation='horizontal', \
ticklocation='bottom', extend='min')
cbar.set_label(r'Inter-model minimum score of the Kernel Ridge Regression')
#---------------------
# Save figure
#---------------------
fignam = 'maps-summarize-score-KR.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# 12s
2023-07-31 17:21:13.363004 Maps summarizing the score of KR: fig File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220602-score-analysis/data2plot-maps-summarize-score-KR-360x180.pckl Figure saved: maps-summarize-score-KR.png CPU times: user 19.8 s, sys: 240 ms, total: 20.1 s Wall time: 19.9 s
%%time
print(datetime.datetime.now())
print('Histogram and interquartile of the inter-model min score')
zfact = 1e6*3600*24 # kmol/m2/s -> mmol/m2/d
#-----------------
# Load data2plot
#-----------------
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
#-----------------
# FIGURE PARAM
#-----------------
nrow = 2
ncol = 2
fsize = (18*cm2in, 16*cm2in) # width height
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# KEYWORDS DICT
#-----------------
kwhist = dict(bins=50, range=[0, 1], cumulative=True, density=True, color='black')
kwlines = dict(color='gray', ls='--')
#---------------------
# Plot histograms
#---------------------
for isimu, simu in enumerate(simu_list):
zax = ax.flat[isimu]
zwdata = data2plot['maps'][simu]['min score'].flatten()
pcm = zax.hist(zwdata, **kwhist)
quantile50 = np.nanquantile(zwdata, .5)
quantile20 = np.nanquantile(zwdata, .2)
quantile05 = np.nanquantile(zwdata, .05)
zax.axhline(0.05, **kwlines)
zax.axhline(0.2, **kwlines)
zax.axhline(0.5, **kwlines)
zax.axvline(quantile50, **kwlines)
zax.axvline(quantile20, **kwlines)
zax.axvline(quantile05, **kwlines)
zax.set_title(subnum.pop()+") "+simu.upper()+"\n95pct above: %.2f\n80pct above: %.2f\n50pct above: %.2f" %(quantile05, quantile20, quantile50), loc='left')
#
fig.suptitle('Histogram and interquartile of the inter-model min score', y=1.03)
fig.tight_layout()
# ca. 3s
2022-12-05 17:58:22.602333 Histogram and interquartile of the inter-model min score File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220602-score-analysis/data2plot-maps-summarize-score-KR-360x180.pckl CPU times: user 2.99 s, sys: 653 ms, total: 3.64 s Wall time: 1.89 s
%%time
print(datetime.datetime.now())
print('Maps models with lowest score')
#-----------------
# Load data2plot
#-----------------
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
zfact = 1e6*3600*24 # kmol/m2/s -> mmol/m2/d
#-----------------
# FIGURE PARAM
#-----------------
nrow = 2
ncol = 2
fsize = (18*cm2in, 16*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# KEYWORDS DICT
#-----------------
#_________________
# Create new cmap
cmap = cm.get_cmap('rainbow', 256)
newcolors = cmap(np.linspace(0, 1, 5))
cmap = ListedColormap(newcolors[:])
cmap.set_under(color='silver')
cmap.set_over(color='silver')
cmap.set_bad(color='silver', alpha=0)
kwmap = {'vmin':-.5, 'vmax':4.5, 'cmap':cmap, \
'transform':ccrs.PlateCarree() }
#---------------------
# Plot score maps
#---------------------
for isimu, simu in enumerate(simu_list):
zax = ax.flat[isimu]
X = data2plot['maps'][simu]['X']
Y = data2plot['maps'][simu]['Y']
Z = data2plot['maps'][simu]['worse model']
pcm = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
zax.set_title(subnum.pop() + ') ' + simu_name[simu].upper())
#
fig.tight_layout()
#---------------------
# Repositionning axes
#---------------------
for icol, zax in enumerate(ax[1, :]):
zw1 = zax.get_position()
ny0 = zw1.y0 + .4*zw1.height
zax.set_position([zw1.x0, ny0, zw1.width, zw1.height])
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, 1].get_position()
nx0 = zw1.x0 + 0.3*zw1.width
ny0 = zw1.y1 + .3*zw1.height
nw = zw2.x1 - 0.3*zw2.width - nx0
nh = 0.05*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm, cax=cax, orientation='horizontal', \
ticklocation='top', extend='neither', ticks=np.arange(0,5))
cbar.ax.set_xticklabels(data2plot['maps']['picontrol']['model list ordered'])
cbar.set_label(r'Model with minimum score of the Kernel Regression', size='x-large')
#---------------------
# Save figure
#---------------------
fignam = 'maps-models-lowest-score-KR.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# 12s
2022-11-21 17:50:16.020120 Maps models with lowest score File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220602-score-analysis/data2plot-maps-summarize-score-KR-360x180.pckl Figure saved: maps-models-lowest-score-KR.png CPU times: user 6.79 s, sys: 839 ms, total: 7.63 s Wall time: 5.24 s
Inspired by the code in dimred-220826-cluster-analysis.ipynb Reorganize clusters with decreasing fiALK/fiDIC ratio + stippling for area with sum of fiALK and fiDIC lower than 0.6
data2plot_dir = dirroot + 'BIOME/V6/'
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
feature_list = ['siconc', 'talk', 'dissic', 'tos', 'sfcWind', 'sos', 'po4', 'atmco2']
nboldclusters = 10
grid = '360x180'
%%time
print(datetime.datetime.now())
print('Maps clusters: clusters grouping and naming')
# group 1: sea ice
# group 2: DIC/ALK (nb: this relate to the buffer factor Samiento and Gruber, 2006, p. 333)
# group 3: Specific cluster
clusters_old_groups_dicts = {
'NorESM2-LM': {
'group 1': [10],
'group 2': [1, 2, 3, 4, 5, 6, 7, 8, 9]
},
'MPI-ESM1-2-LR': {
'group 1': [10],
'group 2': [1, 2, 3, 4, 5, 6, 7, 8, 9]
},
'CESM2': {
'group 1': [10],
'group 2': [1, 2, 3, 4, 5, 6, 7, 8, 9]
},
'ACCESS-ESM1-5': {
'group 1': [10],
'group 2': [1, 2, 3, 4, 5, 6, 7, 8, 9]
},
'IPSL-CM6A-LR': {
'group 1': [10],
'group 2': [1, 2, 3, 4, 5, 6, 7, 8, 9]
}
}
#-----------------
# CREATE CUSTOM CMAP
#-----------------
cmap_dict = {}
for kkk in clusters_old_groups_dicts.keys():
zw = [cm.get_cmap('cool' , 256)(np.linspace(0.1, 0.3, len(clusters_old_groups_dicts[kkk]['group 1']) )),
cm.get_cmap('plasma_r', 256)(np.linspace(0, .7, len(clusters_old_groups_dicts[kkk]['group 2'])))]
newcolors=np.concatenate(zw)
cmap = ListedColormap(newcolors)
cmap.set_under(color='silver')
cmap.set_over(color='silver')
cmap.set_bad(color='silver', alpha=0)
cmap_dict[kkk] = cmap
#
2023-07-27 20:23:43.773971 Maps clusters: clusters grouping and naming CPU times: user 6.06 ms, sys: 807 µs, total: 6.87 ms Wall time: 5.81 ms
def read_biome_datas(zwfname, zwresolution):
if zwresolution == 1: zwgridshape = (360, 180)
elif zwresolution == 2: zwgridshape = (180, 90)
elif zwresolution == 4: zwgridshape = (90, 45)
else: exit('zwresolution not correct')
zwdata = np.loadtxt(zwfname, dtype=str)
# zwdata = np.genfromtxt(zwfname, dtype=str)
Z = np.array(np.where(zwdata.flatten()=='NA', np.nan, zwdata.flatten()).reshape(zwgridshape), dtype=float).T
# X and Y are adapted for plt.pcolormesh with shading='flat'
zwX = np.arange(-180, 180+zwresolution/2., zwresolution)
zwY = np.arange(-90, 90+zwresolution/2., zwresolution)
X, Y = np.meshgrid(zwX, zwY)
return Z, Y, X
#
%%time
print(datetime.datetime.now())
print('Maps clusters: load and prepare data2plot')
if grid == '360x180': resolution = 1
elif grid == '180x90' : resolution = 2
elif grid == '90x45' : resolution = 4
else: exit('grid not correct')
###################
# Load clusters
###################
clusters_old = {}
for model in model_list:
clusters_old[model] = {}
for simu in simu_list:
fname_cluster = glob.glob(data2plot_dir + 'RPART9/' + model + '_' + simu + \
'_grid' + grid + '_K_*_BIOME_RPART_9_opt_3D_SI.txt')[0]
clusters_old[model][simu] = {'X': read_biome_datas(fname_cluster, resolution)[2],
'Y': read_biome_datas(fname_cluster, resolution)[1],
'Z': read_biome_datas(fname_cluster, resolution)[0]}
print("File loaded: "+fname_cluster)
#
#
###################
# Old clusters composition
###################
print('> Old clusters composition')
clusters_old_composition = {}
for model in model_list:
feature_importances = {}
#------------------
# Read feature importances
#------------------
for simu in simu_list:
fname_fi = dirshared + 'fi-std-' + model + '-' + simu + '-' + grid + '.pckl'
with open(fname_fi, 'rb') as f1: feature_importances[simu] = pickle.load(f1)
print("File loaded: "+fname_fi)
#
#------------------
# Compute cluster compositions
#------------------
clusters_old_composition[model] = {}
for feature in feature_list:
clusters_old_composition[model][feature] = {'avg':[], 'std':[]}
for cluster_num in np.arange(1, nboldclusters+.5):
zw= []
for simu in simu_list:
zwfeat_avg = feature_importances[simu][feature]['avg'].T
zw.extend(zwfeat_avg[np.where( clusters_old[model][simu]['Z']==cluster_num )])
#
clusters_old_composition[model][feature]['avg'].append(np.nanmean(zw))
clusters_old_composition[model][feature]['std'].append(np.nanstd (zw))
#
#
#
###################
###################
###################
# Organizing new clusters in groups
# group 1: sea ice importance decreasing
# group 2: ratio DIC importance to ALK importance decreasing
###################
###################
###################
print('> Organizing new clusters in groups')
clusters_new1_groups_dicts = {}
for model in model_list:
clusters_new1_groups_dicts[model] = {}
# group 1
zwclusters = clusters_old_groups_dicts[model]['group 1']
zwsiconc = clusters_old_composition[model]['siconc']['avg']
# remove 1e20 from the index list
zwidx = np.array(np.array(zwclusters)[np.array(zwclusters)<1e20], dtype=int)
# select zwsiconc indexes
zwsiconc2 = []
for iii in zwidx: zwsiconc2.append(zwsiconc[iii-1])
# sort
zwidx_sorted=[iii for _,iii in sorted(zip(zwsiconc2,zwidx))][::-1]
# save
clusters_new1_groups_dicts[model]['group 1'] = zwidx_sorted
# group 2
zwclusters = clusters_old_groups_dicts[model]['group 2']
zwdic = clusters_old_composition[model]['dissic']['avg']
zwalk = clusters_old_composition[model]['talk']['avg']
# compute dic / alk
zwrat=[]
for iii in np.arange(len(zwdic)): zwrat.append(zwdic[iii]/zwalk[iii])
# remove 1e20 from the index list
zwidx = np.array(np.array(zwclusters)[np.array(zwclusters)<1e20], dtype=int)
# select zwrat indexes
zwrat2 = []
for iii in zwidx: zwrat2.append(zwrat[iii-1])
# sort
zwidx_sorted=[iii for _,iii in sorted(zip(zwrat2,zwidx))]
# save
clusters_new1_groups_dicts[model]['group 2'] = zwidx_sorted
#
###################
# Create new to old clusters num dicts and vice versa
###################
print('> Create new1 to old clusters num dicts')
clusters_new1_to_old_dicts = {}
for key1, item1 in clusters_new1_groups_dicts.items():
clusters_new1_to_old_dicts[key1] = {}
icnt = 1
for key2, item2 in item1.items():
for ccc in item2:
clusters_new1_to_old_dicts[key1][str(icnt)]=ccc
icnt+=1
#
#
#
clusters_old_to_new1_dicts = {}
for key1, value1 in clusters_new1_to_old_dicts.items():
clusters_old_to_new1_dicts[key1] = {}
for key2, value2 in value1.items():
clusters_old_to_new1_dicts[key1][str(value2)] = int(key2)
#
#
nbnew1clusters_dict = {}
for model in model_list: nbnew1clusters_dict[model] = len(clusters_new1_to_old_dicts[model].keys())
###################
# New clusters
###################
print('> New clusters')
clusters_new1 = {}
for model in model_list:
clusters_new1[model] = {}
for simu in simu_list:
clusters_new1[model][simu] = {}
#--------------------
# load X, Y, Z and dict
#--------------------
zwX = clusters_old[model][simu]['X']
zwY = clusters_old[model][simu]['Y']
zwZ = clusters_old[model][simu]['Z']
land = np.isnan(zwZ)
zw_new_to_old = clusters_new1_to_old_dicts[model]
#--------------------
# Compute new clusters
#--------------------
nwZ = np.zeros_like(zwZ)
for zwnew, zwold in zw_new_to_old.items():
nwZ = np.where(zwZ == zwold, float(zwnew), nwZ)
#
nwZ[land]=np.nan
#--------------------
# Save
#--------------------
clusters_new1[model][simu]['X'] = zwX
clusters_new1[model][simu]['Y'] = zwY
clusters_new1[model][simu]['Z'] = nwZ
#
#
###################
# New clusters composition
###################
print('> New clusters composition')
clusters_new1_composition = {}
for model in model_list:
feature_importances = {}
clusters_new1_composition[model] = {}
#------------------
# Read feature importances
#------------------
for simu in simu_list:
fname_fi = dirshared + 'fi-std-' + model + '-' + simu + '-' + grid + '.pckl'
with open(fname_fi, 'rb') as f1: feature_importances[simu] = pickle.load(f1)
print("File loaded: "+fname_fi)
#
#------------------
# Compute cluster compositions
#------------------
clusters_new1_composition[model] = {}
for feature in feature_list:
clusters_new1_composition[model][feature] = {'avg':[], 'std':[]}
for cluster_num in np.arange(1, nbnew1clusters_dict[model]+.5):
zw= []
for simu in simu_list:
zwfeat_avg = feature_importances[simu][feature]['avg'].T
zw.extend(zwfeat_avg[np.where( clusters_new1[model][simu]['Z']==cluster_num )])
#
clusters_new1_composition[model][feature]['avg'].append(np.nanmean(zw))
clusters_new1_composition[model][feature]['std'].append(np.nanstd (zw))
#
#
#
###################
# New clusters ocean area
###################
print("> New clusters ocean area")
clusters_new1_ocean_area = {}
earthradius = 6371 # km
for model in model_list:
clusters_new1_ocean_area[model] = {}
for simu in simu_list:
clusters_new1_ocean_area[model][simu] = []
zwX = clusters_new1[model][simu]['X']
zwY = clusters_new1[model][simu]['Y']
zwZ = clusters_new1[model][simu]['Z']
zwcellarea = (np.sin( np.deg2rad(zwY[1:, :-1]) ) - np.sin(np.deg2rad(zwY[:-1, :-1]))) * \
(np.deg2rad(zwX[:-1, 1:]) - np.deg2rad(zwX[:-1, :-1])) * earthradius**2
totoceanarea = np.nansum( np.where(np.isfinite(zwZ), zwcellarea, np.nan) )
for cluster_num in np.arange(1, nbnew1clusters_dict[model]+.5):
clusters_new1_ocean_area[model][simu].append( np.nansum( np.where(zwZ == cluster_num, zwcellarea, np.nan) ) / totoceanarea)
#
#
#
###################
###################
###################
# Organizing new clusters in groups
# and Create new to old clusters num dicts
# group 1: 0
# group 2: 1 for cluster with sum of DIC and ALK importance below treshold
###################
###################
###################
treshold = 0.6
print('> Organizing new clusters in groups and create new to old clusters dicts')
clusters_new2_to_old_dicts = {}
for model in model_list:
clusters_new2_to_old_dicts[model] = {'0':[], '1':[]}
# group 1
clusters_new2_to_old_dicts[model]['0'].extend( clusters_old_groups_dicts[model]['group 1'] )
# group 2
zwclusters = clusters_old_groups_dicts[model]['group 2']
zwdic = clusters_old_composition[model]['dissic']['avg']
zwalk = clusters_old_composition[model]['talk']['avg']
# compute dic + alk
zwsum=[]
for iii in np.arange(len(zwdic)): zwsum.append(zwdic[iii] + zwalk[iii])
# remove 1e20 from the index list and select clusters iidx
zwidx = np.array(np.array(zwclusters)[np.array(zwclusters)<1e20], dtype=int)
# sort clusters
for iii in zwidx:
if zwsum[iii-1] <= treshold: clusters_new2_to_old_dicts[model]['1'].append(iii)
elif zwsum[iii-1] > treshold: clusters_new2_to_old_dicts[model]['0'].append(iii)
else: exit('something weird ini zwsum')
#
#
nbnew2clusters_dict = {}
for model in model_list: nbnew2clusters_dict[model] = len(clusters_new2_to_old_dicts[model].keys())
###################
# Create new2 to new1 clusters num dicts
###################
clusters_new2_to_new1_dicts = {}
for kmodelname, vclusters_new2_to_old in clusters_new2_to_old_dicts.items():
clusters_new2_to_new1_dicts[kmodelname] = {}
for knew2, vold in vclusters_new2_to_old.items():
clusters_new2_to_new1_dicts[kmodelname][knew2] = []
for zwold in vold:
clusters_new2_to_new1_dicts[kmodelname][knew2].append(clusters_old_to_new1_dicts[kmodelname][str(zwold)])
#
#
#
###################
# New clusters
###################
print('> New clusters')
clusters_new2 = {}
for model in model_list:
clusters_new2[model] = {}
for simu in simu_list:
clusters_new2[model][simu] = {}
#--------------------
# load X, Y, Z and dict
#--------------------
zwX = clusters_old[model][simu]['X']
zwY = clusters_old[model][simu]['Y']
zwZ = clusters_old[model][simu]['Z']
land = np.isnan(zwZ)
zw_new_to_old = clusters_new2_to_old_dicts[model]
#--------------------
# Compute new clusters
#--------------------
nwZ = np.zeros_like(zwZ)
for zwnew, zwold_list in zw_new_to_old.items():
for zwold in zwold_list:
nwZ = np.where(zwZ == zwold, float(zwnew), nwZ)
#
#
nwZ[land]=np.nan
#--------------------
# Save
#--------------------
clusters_new2[model][simu]['X'] = .5*(zwX[:, 1:]+zwX[:, :-1])[:-1, :]
clusters_new2[model][simu]['Y'] = .5*(zwY[1:, :]+zwY[:-1, :])[:, :-1]
clusters_new2[model][simu]['Z'] = nwZ
#
#
# ca. 8s
2023-07-27 20:23:55.920979 Maps clusters: load and prepare data2plot File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/ACCESS-ESM1-5_picontrol_grid360x180_K_80_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/ACCESS-ESM1-5_historical_grid360x180_K_80_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/ACCESS-ESM1-5_ssp126_grid360x180_K_80_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/ACCESS-ESM1-5_ssp585_grid360x180_K_80_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/CESM2_picontrol_grid360x180_K_105_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/CESM2_historical_grid360x180_K_105_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/CESM2_ssp126_grid360x180_K_105_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/CESM2_ssp585_grid360x180_K_105_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/IPSL-CM6A-LR_picontrol_grid360x180_K_85_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/IPSL-CM6A-LR_historical_grid360x180_K_85_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/IPSL-CM6A-LR_ssp126_grid360x180_K_85_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/IPSL-CM6A-LR_ssp585_grid360x180_K_85_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/MPI-ESM1-2-LR_picontrol_grid360x180_K_93_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/MPI-ESM1-2-LR_historical_grid360x180_K_93_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/MPI-ESM1-2-LR_ssp126_grid360x180_K_93_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/MPI-ESM1-2-LR_ssp585_grid360x180_K_93_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/NorESM2-LM_picontrol_grid360x180_K_91_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/NorESM2-LM_historical_grid360x180_K_91_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/NorESM2-LM_ssp126_grid360x180_K_91_BIOME_RPART_9_opt_3D_SI.txt File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/V6/RPART9/NorESM2-LM_ssp585_grid360x180_K_91_BIOME_RPART_9_opt_3D_SI.txt > Old clusters composition File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp585-360x180.pckl > Organizing new clusters in groups > Create new1 to old clusters num dicts > New clusters > New clusters composition File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp585-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-picontrol-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-historical-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp126-360x180.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp585-360x180.pckl > New clusters ocean area > Organizing new clusters in groups and create new to old clusters dicts > New clusters CPU times: user 5.32 s, sys: 65.3 ms, total: 5.38 s Wall time: 6.43 s
%%time
print(datetime.datetime.now())
print('Maps clusters: plot')
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (20*cm2in, 18*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
biomes_fay_file='/mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/Time_Varying_Biomes.nc'
zwds = xr.open_dataset(biomes_fay_file)
biomes_fay = {'X': zwds['MeanBiomes'].T['lon'].values,
'Y': zwds['MeanBiomes'].T['lat'].values,
'Z': zwds['MeanBiomes'].T.values}
nbiomes_fay = np.nanmax(biomes_fay['Z'])
##################
##################
#
# MAPS
#
##################
##################
#---------------------
# Plot
#---------------------
irow = 0
pcm = np.zeros_like(ax)
for model in model_list:
icol = 0
kwmap = dict(cmap=cmap_dict[model], vmin=.5, \
vmax=nbnew1clusters_dict[model]+.5, \
transform=ccrs.PlateCarree() )
dens = 4
kwmap2 = dict(transform=ccrs.PlateCarree(), colors='none', hatches=[dens*'.', dens*'.'])
for simu in simu_list:
zax = ax[irow, icol]
X = clusters_new1[model][simu]['X']
Y = clusters_new1[model][simu]['Y']
Z = clusters_new1[model][simu]['Z']
zwpcm = zax.pcolormesh(X, Y, Z, **kwmap)
if simu == 'historical':
zax.contour(biomes_fay['X'], biomes_fay['Y'], biomes_fay['Z'], \
levels=np.arange(1, nbiomes_fay+.5), colors='k', \
transform=ccrs.PlateCarree(), linewidths=1.)
X2 = clusters_new2[model][simu]['X']
Y2 = clusters_new2[model][simu]['Y']
Z2 = clusters_new2[model][simu]['Z']
zax.contourf(X2, Y2, Z2, levels = [0.5, 1.5], **kwmap2)
zax.set_global()
#zax.set_title('%s, %s' %(model.upper(), simu.upper()))
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
pcm[irow, icol] = zwpcm
icol+=1
#
irow+=1
#
fig.tight_layout()
#---------------------
# Move axis
#---------------------
dy = .2*ax[0, 0].get_position().height
for iaxrow, axrow in enumerate(ax[1:]):
for zax in axrow:
zw = zax.get_position()
ny0 = zw.y0 + (iaxrow+1)*dy
zax.set_position([zw.x0, ny0, zw.width, zw.height])
#
#
##################
##################
#
# CLUSTER BARPLOT
#
##################
##################
#-----------------
# KEYWORDS DICT
#-----------------
nbfeature = len(feature_list)
www = 1/(nbfeature+3)
kwbar = dict(align='edge', width=www)
featurecolors = {
'siconc' :'tab:cyan',
'talk' :'tab:orange',
'dissic' :'tab:blue',
'tos' :'tab:green',
'sfcWind':'tab:red',
'sos' :'tab:brown',
'po4' :'tab:olive',
'atmco2' :'tab:gray'
}
simucolors = {
'picontrol':'black',
'historical':'darkgray',
'ssp126':'dodgerblue',
'ssp585':'orangered'
}
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
#---------------------
# Plot
#---------------------
axbarplot1 = []
axbarplot2 = []
cbar = []
for imodel, model in enumerate(model_list):
#_________________
# relative importances
zw = ax[imodel, -1].get_position()
nx0 = zw.x1 + 0.3*zw.width
ny0 = zw.y0+0.2*zw.height
nh = 0.8*zw.height
nw = zw.width
zax1 = fig.add_axes([nx0, ny0, nw, nh])
axbarplot1.append(zax1)
www=0.7
xxx=np.arange(1, nbnew1clusters_dict[model]+1)
www1=www/nbfeature
xxx1=xxx-0.35+www1/2.
bbb=np.zeros_like(xxx)
for feature in feature_list:
hhh = clusters_new1_composition[model][feature]['avg']
hhh = np.where(np.isnan(hhh), 0, hhh)
zax1.bar(xxx, hhh, bottom=bbb, label=feature_name[feature],
#yerr=cluster_composition[model][feature]['std'],
width=www, align='center', color=featurecolors[feature])
zax1.errorbar(xxx1, bbb+np.array(hhh)/2,
yerr=clusters_new1_composition[model][feature]['std'],
fmt='none', ecolor='k', elinewidth=1., capsize=2)
xxx1=xxx1+www1
bbb = bbb+hhh
#
hdls1 = zax1.get_legend_handles_labels()
#_________________
# add colorbar relative importances
zpcm = pcm[imodel, 0]
zw = zax1.get_position()
nx0 = zw.x0 + .015*zw.width
nw = zw.x1 - .015*zw.width - nx0
ny0 = zw.y0 + .05*zw.height
nh = .05 * nw
zwcax = fig.add_axes([nx0, ny0, nw, nh])
zwcbar = plt.colorbar(zpcm, cax=zwcax, orientation='horizontal', \
extend='neither', \
ticks=[])
cbar.append(zwcbar)
#_________________
# add stippling on colorbar relative importances
zw = zwcax.get_position()
ny0 = zw.y0
nh = zw.height
for zwcluster in clusters_new2_to_new1_dicts[model]['1']:
nw = 1/nbnew1clusters_dict[model] * zw.width
nx0 = zw.x0 + (zwcluster-1)*nw
zaxstippling=fig.add_axes([nx0, ny0, nw, nh])
ppp=mpl.patches.Rectangle((0, 0), 1, 1, hatch=dens*'.', fill=False, lw=0)
zaxstippling.add_patch(ppp)
zaxstippling.set_axis_off()
#
#_________________
# set x and y axis relative importances
xtickname = []
zax1.set_xticks(xxx)
# if imodel==len(model_list)-1:
# for iii in xxx: xtickname.append('cluster '+str(iii))
# zax1.set_xticklabels(xtickname, rotation=40, horizontalalignment='right')
# else: zax1.set_xticklabels('')
for iii in xxx: xtickname.append(str(iii))
zax1.set_xticklabels(xtickname)
zax1.set_ylim((-.3, 1.1))
zax1.set_xlim((xxx[0]-www, xxx[-1]+www))
#zax1.set_title(model.upper()+', relative importances')
zax1.axhline(treshold, lw=1, c='k', ls='--')
#_________________
# Ocean areas
zw = zax1.get_position()
nx0 = zw.x1 + 0.1*zw.width
ny0 = zw.y0
nh = zw.height
nw = zw.width
zax2 = fig.add_axes([nx0, ny0, nw, nh])
axbarplot2.append(zax2)
xxx=np.arange(1, nbnew1clusters_dict[model]+1)
www=0.7/3
for simu in simu_list:
zax2.bar(xxx, clusters_new1_ocean_area[model][simu], label=simu_name[simu],
width=www, align='edge', color=simucolors[simu])
xxx=xxx+www
#
zax2.yaxis.tick_right()
hdls2 = zax2.get_legend_handles_labels()
#_________________
# add colorbar Ocean areas
zpcm = pcm[imodel, 0]
zw = zax2.get_position()
nx0 = zw.x0 + .015*zw.width
nw = zw.x1 - .015*zw.width - nx0
ny0 = zw.y0 + .05*zw.height
nh = .05 * nw
zwcax = fig.add_axes([nx0, ny0, nw, nh])
zwcbar = plt.colorbar(zpcm, cax=zwcax, orientation='horizontal', \
extend='neither', \
ticks=[])
cbar.append(zwcbar)
#_________________
# add stippling on colorbar Ocean areas
zw = zwcax.get_position()
ny0 = zw.y0
nh = zw.height
for zwcluster in clusters_new2_to_new1_dicts[model]['1']:
nw = 1/nbnew1clusters_dict[model] * zw.width
nx0 = zw.x0 + (zwcluster-1)*nw
zaxstippling=fig.add_axes([nx0, ny0, nw, nh])
ppp=mpl.patches.Rectangle((0, 0), 1, 1, hatch=dens*'.', fill=False, lw=0)
zaxstippling.add_patch(ppp)
zaxstippling.set_axis_off()
#
#_________________
# set x and y axis Ocean areas
xxx0 = np.arange(1, nbnew1clusters_dict[model]+1)
xtickpos = .5*(xxx-xxx0)+xxx0
xtickname = []
num=0
zax2.set_xticks(xtickpos)
# if imodel==len(model_list)-1:
# for iii in xtickpos:
# xtickname.append('cluster '+str(num+1))
# num+=1
# #
# zax2.set_xticklabels(xtickname, rotation = 40, horizontalalignment='right')
# else: zax2.set_xticklabels('')
for iii in xtickpos:
xtickname.append(str(num+1))
num+=1
#
zax2.set_xticklabels(xtickname)
zax2.set_ylim((-.164, .6))
zax2.set_xlim((xtickpos[0]-3*www, xtickpos[-1]+3*www))
#zax2.set_title(model.upper()+', ocean areas')
#
# remove cbar borders
for zwcbar in cbar: zwcbar.outline.set_visible(False)
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') ' + zax.get_title()
zax.set_title(zwtitle, loc='left')
#
zwtitle = subnum.pop() + ') '+ axbarplot1[irow].get_title()
axbarplot1[irow].set_title(zwtitle, loc='left')
zwtitle = subnum.pop() + ') '+ axbarplot2[irow].get_title()
axbarplot2[irow].set_title(zwtitle, loc='left')
#
#---------------------
# Add legends
#---------------------
kwleg = dict(handlelength=1, handleheight=.7, handletextpad=.5, columnspacing=1., \
loc='lower center')
zw = axbarplot1[0].get_position()
nx0 = zw.x1 - 0.5*zw.width
ny0 = zw.y1+0.3*zw.height
zax = fig.add_axes([nx0, ny0, 0, 0])
zax.axis('off')
leg=zax.legend(handles=hdls1[0], labels=hdls1[1], ncol=2, **kwleg)
leg.set_title('Cluster compositions: ', prop=dict(weight='bold'))
zw = axbarplot2[0].get_position()
nx0 = zw.x1 - 0.5*zw.width
ny0 = zw.y1+0.3*zw.height
zax = fig.add_axes([nx0, ny0, 0, 0])
zax.axis('off')
leg=zax.legend(handles=hdls2[0], labels=hdls2[1], ncol=2, **kwleg)
leg.set_title('Cluster ocean areas ', prop=dict(weight='bold'))
#---------------------
# Add text to describe cluster
#---------------------
zw = axbarplot1[-1].get_position()
ny0 = zw.y0 - .4*zw.height
nx0 = zw.x0 + 1/11.*zw.width
fig.text(nx0, ny0, 'Sea-ice\ncluster', ha='right', va='top')
nx0 = zw.x0 + 1.5/11.*zw.width
fig.text(nx0, ny0, 'ALK\nclusters', ha='left', va='top')
nx0 = zw.x0 + 10.5/11*zw.width
fig.text(nx0, ny0, 'DIC\nclusters', ha='right', va='top')
zw = axbarplot2[-1].get_position()
ny0 = zw.y0 - .4*zw.height
nx0 = zw.x0 + 1/11.*zw.width
fig.text(nx0, ny0, 'Sea-ice\ncluster', ha='right', va='top')
nx0 = zw.x0 + 1.5/11.*zw.width
fig.text(nx0, ny0, 'ALK\nclusters', ha='left', va='top')
nx0 = zw.x0 + 10.5/11*zw.width
fig.text(nx0, ny0, 'DIC\nclusters', ha='right', va='top')
#---------------------
# Add suptitle
#---------------------
zw = ax[0, 0].get_position()
nx0 = zw.x0
ny0 = zw.y1+0.8*zw.height
fig.suptitle('Clusters organized with decreasing fiALK/fiDIC ratio', y=ny0, x=nx0, \
size='large', ha='left', va='bottom')
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .1*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .2*zw1.height
fig.text(nx0, ny0, vsimu.upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-clusters.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 50s
2022-11-21 17:21:36.278445 Maps clusters: plot Figure saved: maps-clusters.png CPU times: user 57.1 s, sys: 5.3 s, total: 1min 2s Wall time: 48.2 s
%%time
print(datetime.datetime.now())
print('Maps clusters: plot summary')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['historical', 'ssp126', 'ssp585']
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (14*cm2in, 17*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
biomes_fay_file='/mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/Time_Varying_Biomes.nc'
zwds = xr.open_dataset(biomes_fay_file)
biomes_fay = {'X': zwds['MeanBiomes'].T['lon'].values,
'Y': zwds['MeanBiomes'].T['lat'].values,
'Z': zwds['MeanBiomes'].T.values}
nbiomes_fay = np.nanmax(biomes_fay['Z'])
##################
##################
#
# MAPS
#
##################
##################
#---------------------
# Plot
#---------------------
irow = 0
pcm = np.zeros_like(ax)
for model in model_list:
icol = 0
kwmap = dict(cmap=cmap_dict[model], vmin=.5, \
vmax=nbnew1clusters_dict[model]+.5, \
transform=ccrs.PlateCarree() )
dens = 4
kwmap2 = dict(transform=ccrs.PlateCarree(), colors='none', hatches=[dens*'.', dens*'.'])
for simu in simu_list:
zax = ax[irow, icol]
X = clusters_new1[model][simu]['X']
Y = clusters_new1[model][simu]['Y']
Z = clusters_new1[model][simu]['Z']
zwpcm = zax.pcolormesh(X, Y, Z, **kwmap)
if simu == 'historical':
zax.contour(biomes_fay['X'], biomes_fay['Y'], biomes_fay['Z'], \
levels=np.arange(1, nbiomes_fay+.5), colors='k', \
transform=ccrs.PlateCarree(), linewidths=1.)
X2 = clusters_new2[model][simu]['X']
Y2 = clusters_new2[model][simu]['Y']
Z2 = clusters_new2[model][simu]['Z']
zax.contourf(X2, Y2, Z2, levels = [0.5, 1.5], **kwmap2)
zax.set_global()
#zax.set_title('%s, %s' %(model.upper(), simu.upper()))
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
pcm[irow, icol] = zwpcm
icol+=1
#
irow+=1
#
fig.tight_layout()
#---------------------
# Adjust maps axis
#---------------------
for iaxrow, axrow in enumerate(ax):
for zax in axrow:
zw = zax.get_position()
zax.set_position([zw.x0, zw.y0, 1.06*zw.width, 1.07*zw.height])
#
#
dy = .25*ax[0, 0].get_position().height
for iaxrow, axrow in enumerate(ax[1:]):
for zax in axrow:
zw = zax.get_position()
ny0 = zw.y0 + (iaxrow+1)*dy
zax.set_position([zw.x0, ny0, zw.width, zw.height])
#
#
##################
##################
#
# CLUSTER BARPLOT
#
##################
##################
#-----------------
# KEYWORDS DICT
#-----------------
nbfeature = len(feature_list)
www = 1/(nbfeature+3)
kwbar = dict(align='edge', width=www)
featurecolors = {
'siconc' :'tab:cyan',
'talk' :'tab:orange',
'dissic' :'tab:blue',
'tos' :'tab:green',
'sfcWind':'tab:red',
'sos' :'tab:brown',
'po4' :'tab:olive',
'atmco2' :'tab:gray'
}
simucolors = {
'picontrol':'black',
'historical':'darkgray',
'ssp126':'dodgerblue',
'ssp585':'orangered'
}
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
#---------------------
# Plot
#---------------------
axbarplot1 = []
cbar = []
for imodel, model in enumerate(model_list):
#_________________
# relative importances
zw = ax[imodel, -1].get_position()
nx0 = zw.x1 + 0.15*zw.width
ny0 = zw.y0+0.2*zw.height
nh = 0.8*zw.height
nw = zw.width
zax1 = fig.add_axes([nx0, ny0, nw, nh])
axbarplot1.append(zax1)
www=0.7
xxx=np.arange(1, nbnew1clusters_dict[model]+1)
www1=www/nbfeature
xxx1=xxx-0.35+www1/2.
bbb=np.zeros_like(xxx)
for feature in feature_list:
hhh = clusters_new1_composition[model][feature]['avg']
hhh = np.where(np.isnan(hhh), 0, hhh)
zax1.bar(xxx, hhh, bottom=bbb, label=feature_name[feature],
width=www, align='center', color=featurecolors[feature])
xxx1=xxx1+www1
bbb = bbb+hhh
#
hdls1 = zax1.get_legend_handles_labels()
#_________________
# add colorbar relative importances
zpcm = pcm[imodel, 0]
zw = zax1.get_position()
nx0 = zw.x0 + .015*zw.width
nw = zw.x1 - .015*zw.width - nx0
ny0 = zw.y0 + .05*zw.height
nh = .05 * nw
zwcax = fig.add_axes([nx0, ny0, nw, nh])
zwcbar = plt.colorbar(zpcm, cax=zwcax, orientation='horizontal', \
extend='neither', \
ticks=[])
cbar.append(zwcbar)
#_________________
# add stippling on colorbar relative importances
zw = zwcax.get_position()
ny0 = zw.y0
nh = zw.height
for zwcluster in clusters_new2_to_new1_dicts[model]['1']:
nw = 1/nbnew1clusters_dict[model] * zw.width
nx0 = zw.x0 + (zwcluster-1)*nw
zaxstippling=fig.add_axes([nx0, ny0, nw, nh])
ppp=mpl.patches.Rectangle((0, 0), 1, 1, hatch=dens*'.', fill=False, lw=0)
zaxstippling.add_patch(ppp)
zaxstippling.set_axis_off()
#
#_________________
# set x and y axis relative importances
xtickname = []
zax1.set_xticks(xxx)
# if imodel==len(model_list)-1:
# for iii in xxx: xtickname.append('cluster '+str(iii))
# zax1.set_xticklabels(xtickname, rotation=40, horizontalalignment='right')
# else: zax1.set_xticklabels('')
for iii in xxx: xtickname.append(str(iii))
zax1.set_xticklabels(xtickname)
zax1.yaxis.set_ticks_position('right')
zax1.set_ylim((-.4, 1.1))
zax1.set_xlim((xxx[0]-www, xxx[-1]+www))
#
# remove cbar borders
for zwcbar in cbar: zwcbar.outline.set_visible(False)
##################
##################
#
# DECORATION
#
##################
##################
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') ' + zax.get_title()
#zax.set_title(zwtitle, loc='left')
zax.annotate(zwtitle, (0.16, 0.72),xycoords='axes fraction')
#
zwtitle = subnum.pop() + ') '+ axbarplot1[irow].get_title()
axbarplot1[irow].annotate(zwtitle, (-.09, 0.85), xycoords='axes fraction')
# axbarplot1[irow].set_title(zwtitle, loc='left')
#
#---------------------
# Add legends
#---------------------
kwleg = dict(handlelength=1, handleheight=.7, handletextpad=.3, columnspacing=.7, \
loc='upper left', frameon=False)
zw = axbarplot1[0].get_position()
nx0 = zw.x1 + 0.1*zw.width
ny0 = zw.y1
zax = fig.add_axes([nx0, ny0, 0, 0])
zax.axis('off')
leg=zax.legend(handles=hdls1[0], labels=hdls1[1], ncol=1, **kwleg)
#leg.set_title('Cluster\ncompositions: ', prop=dict(weight='bold'))
#---------------------
# Add text to describe cluster
#---------------------
kwtext = dict(rotation=45., ha='right', va='top')
zw = axbarplot1[-1].get_position()
ny0 = zw.y0 - .4*zw.height
nx0 = zw.x0 + 1/11.*zw.width
fig.text(nx0, ny0, 'Sea-ice\ncluster', ha='right', va='top')
nx0 = zw.x0 + 1.5/11.*zw.width
fig.text(nx0, ny0, 'ALK\nclusters', ha='left', va='top')
nx0 = zw.x0 + 10.5/11*zw.width
fig.text(nx0, ny0, 'DIC\nclusters', ha='right', va='top')
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='right', va='bottom')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = zw1.y1
nx0 = zw1.x0 + 0.15*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .15*zw1.height
fig.text(nx0, ny0, vsimu.upper(), **kwtxt)
#
zw = axbarplot1[0].get_position()
nx0 = zw.x0 + .5*zw.width
fig.text(nx0, ny0, 'Cluster compositions', ha='center', va='bottom')
#---------------------
# Save figure
#---------------------
fignam = 'maps-clusters-summary.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 50s
2023-06-27 17:08:15.552009 Maps clusters: plot summary Figure saved: maps-clusters-summary.png CPU times: user 36.7 s, sys: 527 ms, total: 37.3 s Wall time: 36.9 s
%%time
print(datetime.datetime.now())
print('Maps clusters: plot summary with picontrol')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (17*cm2in, 17*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
biomes_fay_file='/mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/BIOME/Time_Varying_Biomes.nc'
zwds = xr.open_dataset(biomes_fay_file)
biomes_fay = {'X': zwds['MeanBiomes'].T['lon'].values,
'Y': zwds['MeanBiomes'].T['lat'].values,
'Z': zwds['MeanBiomes'].T.values}
nbiomes_fay = np.nanmax(biomes_fay['Z'])
##################
##################
#
# MAPS
#
##################
##################
#---------------------
# Plot
#---------------------
irow = 0
pcm = np.zeros_like(ax)
for model in model_list:
icol = 0
kwmap = dict(cmap=cmap_dict[model], vmin=.5, \
vmax=nbnew1clusters_dict[model]+.5, \
transform=ccrs.PlateCarree() )
dens = 4
kwmap2 = dict(transform=ccrs.PlateCarree(), colors='none', hatches=[dens*'.', dens*'.'])
for simu in simu_list:
zax = ax[irow, icol]
X = clusters_new1[model][simu]['X']
Y = clusters_new1[model][simu]['Y']
Z = clusters_new1[model][simu]['Z']
zwpcm = zax.pcolormesh(X, Y, Z, **kwmap)
if simu == 'historical':
zax.contour(biomes_fay['X'], biomes_fay['Y'], biomes_fay['Z'], \
levels=np.arange(1, nbiomes_fay+.5), colors='k', \
transform=ccrs.PlateCarree(), linewidths=1.)
X2 = clusters_new2[model][simu]['X']
Y2 = clusters_new2[model][simu]['Y']
Z2 = clusters_new2[model][simu]['Z']
zax.contourf(X2, Y2, Z2, levels = [0.5, 1.5], **kwmap2)
zax.set_global()
#zax.set_title('%s, %s' %(model.upper(), simu.upper()))
zax.coastlines()
# gl = zax.gridlines(crs=ccrs.PlateCarree(), linewidth=1)
pcm[irow, icol] = zwpcm
icol+=1
#
irow+=1
#
fig.tight_layout()
#---------------------
# Adjust maps axis
#---------------------
for iaxrow, axrow in enumerate(ax):
for zax in axrow:
zw = zax.get_position()
zax.set_position([zw.x0, zw.y0, 1.06*zw.width, 1.07*zw.height])
#
#
dy = .25*ax[0, 0].get_position().height
for iaxrow, axrow in enumerate(ax[1:]):
for zax in axrow:
zw = zax.get_position()
ny0 = zw.y0 + (iaxrow+1)*dy
zax.set_position([zw.x0, ny0, zw.width, zw.height])
#
#
##################
##################
#
# CLUSTER BARPLOT
#
##################
##################
#-----------------
# KEYWORDS DICT
#-----------------
nbfeature = len(feature_list)
www = 1/(nbfeature+3)
kwbar = dict(align='edge', width=www)
featurecolors = {
'siconc' :'tab:cyan',
'talk' :'tab:orange',
'dissic' :'tab:blue',
'tos' :'tab:green',
'sfcWind':'tab:red',
'sos' :'tab:brown',
'po4' :'tab:olive',
'atmco2' :'tab:gray'
}
simucolors = {
'picontrol':'black',
'historical':'darkgray',
'ssp126':'dodgerblue',
'ssp585':'orangered'
}
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
#---------------------
# Plot
#---------------------
axbarplot1 = []
cbar = []
for imodel, model in enumerate(model_list):
#_________________
# relative importances
zw = ax[imodel, -1].get_position()
nx0 = zw.x1 + 0.15*zw.width
ny0 = zw.y0+0.2*zw.height
nh = 0.8*zw.height
nw = zw.width
zax1 = fig.add_axes([nx0, ny0, nw, nh])
axbarplot1.append(zax1)
www=0.7
xxx=np.arange(1, nbnew1clusters_dict[model]+1)
www1=www/nbfeature
xxx1=xxx-0.35+www1/2.
bbb=np.zeros_like(xxx)
for feature in feature_list:
hhh = clusters_new1_composition[model][feature]['avg']
hhh = np.where(np.isnan(hhh), 0, hhh)
zax1.bar(xxx, hhh, bottom=bbb, label=feature_name[feature],
width=www, align='center', color=featurecolors[feature])
xxx1=xxx1+www1
bbb = bbb+hhh
#
hdls1 = zax1.get_legend_handles_labels()
#_________________
# add colorbar relative importances
zpcm = pcm[imodel, 0]
zw = zax1.get_position()
nx0 = zw.x0 + .015*zw.width
nw = zw.x1 - .015*zw.width - nx0
ny0 = zw.y0 + .05*zw.height
nh = .05 * nw
zwcax = fig.add_axes([nx0, ny0, nw, nh])
zwcbar = plt.colorbar(zpcm, cax=zwcax, orientation='horizontal', \
extend='neither', \
ticks=[])
cbar.append(zwcbar)
#_________________
# add stippling on colorbar relative importances
zw = zwcax.get_position()
ny0 = zw.y0
nh = zw.height
for zwcluster in clusters_new2_to_new1_dicts[model]['1']:
nw = 1/nbnew1clusters_dict[model] * zw.width
nx0 = zw.x0 + (zwcluster-1)*nw
zaxstippling=fig.add_axes([nx0, ny0, nw, nh])
ppp=mpl.patches.Rectangle((0, 0), 1, 1, hatch=dens*'.', fill=False, lw=0)
zaxstippling.add_patch(ppp)
zaxstippling.set_axis_off()
#
#_________________
# set x and y axis relative importances
xtickname = []
zax1.set_xticks(xxx)
# if imodel==len(model_list)-1:
# for iii in xxx: xtickname.append('cluster '+str(iii))
# zax1.set_xticklabels(xtickname, rotation=40, horizontalalignment='right')
# else: zax1.set_xticklabels('')
for iii in xxx: xtickname.append(str(iii))
zax1.set_xticklabels(xtickname)
zax1.yaxis.set_ticks_position('right')
zax1.set_ylim((-.4, 1.1))
zax1.set_xlim((xxx[0]-www, xxx[-1]+www))
#
# remove cbar borders
for zwcbar in cbar: zwcbar.outline.set_visible(False)
##################
##################
#
# DECORATION
#
##################
##################
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') ' + zax.get_title()
#zax.set_title(zwtitle, loc='left')
zax.annotate(zwtitle, (0.16, 0.72),xycoords='axes fraction')
#
zwtitle = subnum.pop() + ') '+ axbarplot1[irow].get_title()
axbarplot1[irow].annotate(zwtitle, (-.09, 0.85), xycoords='axes fraction')
# axbarplot1[irow].set_title(zwtitle, loc='left')
#
#---------------------
# Add legends
#---------------------
kwleg = dict(handlelength=1, handleheight=.7, handletextpad=.3, columnspacing=.7, \
loc='upper left', frameon=False)
zw = axbarplot1[0].get_position()
nx0 = zw.x1 + 0.1*zw.width
ny0 = zw.y1
zax = fig.add_axes([nx0, ny0, 0, 0])
zax.axis('off')
leg=zax.legend(handles=hdls1[0], labels=hdls1[1], ncol=1, **kwleg)
#leg.set_title('Cluster\ncompositions: ', prop=dict(weight='bold'))
#---------------------
# Add text to describe cluster
#---------------------
kwtext = dict(rotation=45., ha='right', va='top')
zw = axbarplot1[-1].get_position()
ny0 = zw.y0 - .4*zw.height
nx0 = zw.x0 + 1/11.*zw.width
fig.text(nx0, ny0, 'Sea-ice\ncluster', ha='right', va='top')
nx0 = zw.x0 + 1.5/11.*zw.width
fig.text(nx0, ny0, 'ALK\nclusters', ha='left', va='top')
nx0 = zw.x0 + 10.5/11*zw.width
fig.text(nx0, ny0, 'DIC\nclusters', ha='right', va='top')
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='right', va='bottom')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = zw1.y1
nx0 = zw1.x0 + 0.15*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .15*zw1.height
fig.text(nx0, ny0, vsimu.upper(), **kwtxt)
#
zw = axbarplot1[0].get_position()
nx0 = zw.x0 + .5*zw.width
fig.text(nx0, ny0, 'Cluster compositions', ha='center', va='bottom')
#---------------------
# Save figure
#---------------------
fignam = 'maps-clusters-summary-with-picontrol.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 50s
2023-07-27 20:26:43.612267 Maps clusters: plot summary with picontrol Figure saved: maps-clusters-summary-with-picontrol.png CPU times: user 47.9 s, sys: 576 ms, total: 48.5 s Wall time: 48.2 s
%%time
print(datetime.datetime.now())
print(' Maps of FI for each models and each simulations: inputs')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
feature_list = ['dissic', 'talk', 'tos', 'sfcWind', 'siconc', 'sos', 'po4']
grid = '360x180'
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2022-11-21 17:24:22.480548 Maps of FI for each models and each simulations: inputs CPU times: user 0 ns, sys: 966 µs, total: 966 µs Wall time: 681 µs
%%time
print(datetime.datetime.now())
print('Maps of FI for each models and each simulations: prepare data2plot')
data2plot = {}
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
data2plot[vmodel] = {}
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
data2plot[vmodel][vsimu] = {}
print('>>> Load datas')
savedfile = dirshared + 'fi-std-' + vmodel + '-' + vsimu + '-' + grid + '.pckl'
with open(savedfile, 'rb') as f1: loadeddatas = pickle.load(f1)
print("File loaded: "+savedfile)
# Load features
print('>>> Get features')
for ifeature, vfeature in enumerate(feature_list):
data2plot[vmodel][vsimu][vfeature] = {'X': loadeddatas['lon'].T,
'Y': loadeddatas['lat'].T,
'Z': loadeddatas[vfeature]['avg'].T}
#
#
#
# ca. 1s
2022-11-21 17:23:05.577562 Maps of FI for each models and each simulations: prepare data2plot > model: ACCESS-ESM1-5 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp585-360x180.pckl >>> Get features > model: CESM2 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp585-360x180.pckl >>> Get features > model: IPSL-CM6A-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp585-360x180.pckl >>> Get features > model: MPI-ESM1-2-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp585-360x180.pckl >>> Get features > model: NorESM2-LM >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp585-360x180.pckl >>> Get features CPU times: user 76 ms, sys: 53.5 ms, total: 130 ms Wall time: 122 ms
%%time
print(datetime.datetime.now())
print('Maps of FI for each models and each simulations: plot')
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
feature_list = ['dissic', 'talk', 'tos', 'sfcWind', 'siconc', 'sos', 'po4']
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(feature_list)
ncol = len(simu_list)
fsize = (18*cm2in, 20*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('RdYlGn_r', 256)
newcolors = zw(np.linspace(0, 1, 10+1))
cmap1 = ListedColormap(newcolors[:-1])
cmap1.set_under(color='silver')
cmap1.set_over (color=newcolors[-1])
cmap1.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':0.7, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
for ifeature, vfeature in enumerate(feature_list):
zax = ax[ifeature, isimu]
X = data2plot[vmodel][vsimu][vfeature]['X']
Y = data2plot[vmodel][vsimu][vfeature]['Y']
Z = data2plot[vmodel][vsimu][vfeature]['Z']
mmm[ifeature, isimu] = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') '
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, -1].get_position()
nx0 = zw1.x0 + .5*zw1.width
ny0 = zw1.y1 + 1*zw1.height
nw = zw2.x1 - .5*zw2.width - nx0
nh = 0.04*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='max',
label='Relative importances in '+vmodel)
cbar.set_ticks(np.linspace(kwmap['vmin'], kwmap['vmax'], mmm[0, 0].cmap.colors.shape[0]+1))
#---------------------
# Add text on the left: feature names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for ifeature, vfeature in enumerate(feature_list):
zw1 = ax[ifeature, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, feature_name[vfeature].upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-fi-'+vmodel+'.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
#
# ca. 2m30s
2022-11-21 17:25:15.713599 Maps of FI for each models and each simulations: plot > model: ACCESS-ESM1-5 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-fi-ACCESS-ESM1-5.png > model: CESM2 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-fi-CESM2.png > model: IPSL-CM6A-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-fi-IPSL-CM6A-LR.png > model: MPI-ESM1-2-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-fi-MPI-ESM1-2-LR.png > model: NorESM2-LM >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-fi-NorESM2-LM.png CPU times: user 3min 1s, sys: 19.8 s, total: 3min 20s Wall time: 2min 27s
row = each model, col = each simu
%%time
print(datetime.datetime.now())
print('Maps of the sum of FI for DIC and ALK: inputs')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
grid = '360x180'
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2022-11-21 17:28:01.090926 Maps of the sum of FI for DIC and ALK: inputs CPU times: user 2.1 ms, sys: 179 µs, total: 2.28 ms Wall time: 520 µs
%%time
print(datetime.datetime.now())
print('Maps of the sum of FI for DIC and ALK: prepare data2plot')
data2plot = {}
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
data2plot[vmodel] = {}
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
data2plot[vmodel][vsimu] = {}
print('>>> Load datas')
savedfile = dirshared + 'fi-std-' + vmodel + '-' + vsimu + '-' + grid + '.pckl'
with open(savedfile, 'rb') as f1: loadeddatas = pickle.load(f1)
print("File loaded: "+savedfile)
# Load features
print('>>> Get features')
data2plot[vmodel][vsimu] = {'X': loadeddatas['lon'].T,
'Y': loadeddatas['lat'].T,
'Z': loadeddatas['dissic']['avg'].T+loadeddatas['talk']['avg'].T}
#
#
#
# ca. 1s
2022-11-21 17:28:01.101968 Maps of the sum of FI for DIC and ALK: prepare data2plot > model: ACCESS-ESM1-5 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp585-360x180.pckl >>> Get features > model: CESM2 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp585-360x180.pckl >>> Get features > model: IPSL-CM6A-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp585-360x180.pckl >>> Get features > model: MPI-ESM1-2-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp585-360x180.pckl >>> Get features > model: NorESM2-LM >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp585-360x180.pckl >>> Get features CPU times: user 90.8 ms, sys: 51.1 ms, total: 142 ms Wall time: 128 ms
%%time
print(datetime.datetime.now())
print('Maps of the sum of FI for DIC and ALK: plot')
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (16*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('RdYlGn_r', 256)
newcolors = zw(np.linspace(0, 1, 10))
cmap1 = ListedColormap(newcolors)
cmap1.set_under(color='silver')
cmap1.set_over (color='silver')
cmap1.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':1., 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
zax = ax[imodel, isimu]
X = data2plot[vmodel][vsimu]['X']
Y = data2plot[vmodel][vsimu]['Y']
Z = data2plot[vmodel][vsimu]['Z']
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') '
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, -1].get_position()
nx0 = zw1.x0 + .5*zw1.width
ny0 = zw1.y1 + 1*zw1.height
nw = zw2.x1 - .5*zw2.width - nx0
nh = 0.04*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='neither',
label='Sum of the relative importances of DIC and ALK')
cbar.set_ticks(np.linspace(kwmap['vmin'], kwmap['vmax'], mmm[0, 0].cmap.colors.shape[0]+1))
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-sum-fi-dic-alk.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 2m30s
2022-11-21 17:28:01.244445 Maps of the sum of FI for DIC and ALK: plot > model: ACCESS-ESM1-5 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-sum-fi-dic-alk.png CPU times: user 25.7 s, sys: 2.77 s, total: 28.4 s Wall time: 21.2 s
%%time
print(datetime.datetime.now())
print('Histogram')
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (18*cm2in, 20*cm2in) # width height
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# KEYWORDS DICT
#-----------------
kwhist = dict(bins=50, range=[0, 1], density=True, color='black')
kwhist2 = dict(bins=50, range=[0, 1], density=True, color='grey')
kwlines = dict(color='gray', ls='--')
#---------------------
# Plot histograms
#---------------------
for imodel, vmodel in enumerate(model_list):
for isimu, vsimu in enumerate(simu_list):
zax = ax[imodel, isimu]
if vsimu == 'historical':
zwdata = data2plot[vmodel][vsimu]['Z'].flatten()
pcm = zax.hist(zwdata, **kwhist2)
else:
zwdata = data2plot[vmodel][vsimu]['Z'].flatten()
pcm = zax.hist(zwdata, **kwhist)
zwdatahist = data2plot[vmodel]['historical']['Z'].flatten()
pcm = zax.hist(zwdatahist, alpha=.5, **kwhist2)
#
quantile20 = np.nanquantile(zwdata, .2)
quantile50 = np.nanquantile(zwdata, .5)
quantile80 = np.nanquantile(zwdata, .8)
zw= np.nonzero(zwdata > .5)
pctabove50 = len(zw[0])/len(np.where(zwdata==zwdata)[0])*100
#zax.axhline(0.2, **kwlines)
#zax.axvline(quantile20, **kwlines)
zax.set_title(subnum.pop()+") %.0f%s > 0.5\n80%s -- 50%s -- 20%s\n%.2f -- %.2f -- %.2f"
%(pctabove50, "%", "%", "%", "%", quantile20, quantile50, quantile80), loc='left')
#
fig.tight_layout()
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .3*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .9*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'histogram-sum-fi-dic-alk.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 7s
2022-11-21 17:28:24.951407 Histogram
/opt/conda/lib/python3.7/site-packages/numpy/lib/histograms.py:839: RuntimeWarning: invalid value encountered in greater_equal keep = (tmp_a >= first_edge) /opt/conda/lib/python3.7/site-packages/numpy/lib/histograms.py:840: RuntimeWarning: invalid value encountered in less_equal keep &= (tmp_a <= last_edge) /opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in greater
Figure saved: histogram-sum-fi-dic-alk.png CPU times: user 6.87 s, sys: 186 ms, total: 7.06 s Wall time: 6.84 s
row = each model, col = each simu
%%time
print(datetime.datetime.now())
print('Maps of the ratio of FIDIC and FIALK: inputs')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['picontrol', 'historical', 'ssp126', 'ssp585']
grid = '360x180'
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2023-06-29 22:06:59.072235 Maps of the ratio of FIDIC and FIALK: inputs CPU times: user 60 µs, sys: 2 µs, total: 62 µs Wall time: 69.9 µs
%%time
print(datetime.datetime.now())
print('Maps of the ratio of FIDIC and FIALK: prepare data2plot')
data2plot = {}
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
data2plot[vmodel] = {}
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
data2plot[vmodel][vsimu] = {}
print('>>> Load datas')
savedfile = dirshared + 'fi-std-' + vmodel + '-' + vsimu + '-' + grid + '.pckl'
with open(savedfile, 'rb') as f1: loadeddatas = pickle.load(f1)
print("File loaded: "+savedfile)
# Load features
print('>>> Get features')
data2plot[vmodel][vsimu] = {'X': loadeddatas['lon'].T,
'Y': loadeddatas['lat'].T,
'Z': loadeddatas['talk']['avg'].T/loadeddatas['dissic']['avg'].T}
#
#
# ca. 1s
2023-06-29 22:07:03.484733 Maps of the ratio of FIDIC and FIALK: prepare data2plot > model: ACCESS-ESM1-5 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-ACCESS-ESM1-5-ssp585-360x180.pckl >>> Get features > model: CESM2 >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-CESM2-ssp585-360x180.pckl >>> Get features > model: IPSL-CM6A-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-IPSL-CM6A-LR-ssp585-360x180.pckl >>> Get features > model: MPI-ESM1-2-LR >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-MPI-ESM1-2-LR-ssp585-360x180.pckl >>> Get features > model: NorESM2-LM >> simu: picontrol >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-picontrol-360x180.pckl >>> Get features >> simu: historical >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-historical-360x180.pckl >>> Get features >> simu: ssp126 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp126-360x180.pckl >>> Get features >> simu: ssp585 >>> Load datas File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/SHARED-DATAS/fi-std-NorESM2-LM-ssp585-360x180.pckl >>> Get features CPU times: user 31.9 ms, sys: 23.9 ms, total: 55.8 ms Wall time: 55.3 ms
<timed exec>:24: RuntimeWarning: divide by zero encountered in true_divide
%%time
print(datetime.datetime.now())
print('Maps of the ratio of FIDIC and FIALK: plot')
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list) + 2
fsize = (24*cm2in, 12*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('RdYlGn_r', 256)
newcolors = zw(np.linspace(0, 1, 10+1))
cmap1 = ListedColormap(newcolors[:-1])
cmap1.set_under(color='silver')
cmap1.set_over (color=newcolors[-1])
cmap1.set_bad(color='silver', alpha=0)
zw = cm.get_cmap('RdBu_r', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap2 = ListedColormap(newcolors[1:-1])
cmap2.set_under(color=newcolors[0])
cmap2.set_over (color=newcolors[-1])
cmap2.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':2., 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2 = {'vmin':-2, 'vmax':2., 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
zax = ax[imodel, isimu]
X = data2plot[vmodel][vsimu]['X']
Y = data2plot[vmodel][vsimu]['Y']
Z = data2plot[vmodel][vsimu]['Z']
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
#
for isimu2, vsimu2 in enumerate(['ssp126', 'ssp585'], start=len(simu_list)):
print('>> simu: ' + vsimu2)
zax = ax[imodel, isimu2]
X = data2plot[vmodel][vsimu2]['X']
Y = data2plot[vmodel][vsimu2]['Y']
Z = data2plot[vmodel][vsimu2]['Z'] - data2plot[vmodel]['historical']['Z']
mmm[imodel, isimu2] = zax.pcolormesh(X, Y, Z, **kwmap2)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ')'
zax.annotate(zwtitle, (0.14, 0.72),xycoords='axes fraction')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 1].get_position()
zw2 = ax[0, len(simu_list)-2].get_position()
nx0 = zw1.x0
ny0 = zw1.y1 + .7*zw1.height
nw = zw2.x1 - nx0
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, 0], cax=cax, ticklocation='top',
orientation='horizontal', extend='max',
label='Ratio of the relative importances of DIC and ALK')
cbar.set_ticks(np.linspace(kwmap['vmin'], kwmap['vmax'], int(mmm[0, 0].cmap.colors.shape[0]/2)+1))
zw1 = ax[0, len(simu_list)].get_position()
zw2 = ax[0, len(simu_list)+1].get_position()
nx0 = zw1.x0
ny0 = zw1.y1 + .7*zw1.height
nw = zw2.x1 - nx0
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='both',
label='Change relative to historical')
cbar.set_ticks(np.linspace(kwmap2['vmin'], kwmap2['vmax'], int(mmm[0, -1].cmap.colors.shape[0]/2)+1))
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='right', va='bottom')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = zw1.y1
nx0 = zw1.x0 + 0.15*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='medium', rotation='horizontal', weight='normal', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .3*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
for isimu2, vsimu2 in enumerate(['ssp126', 'ssp585'], start=len(simu_list)):
zw1 = ax[0, isimu2].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .3*zw1.height
fig.text(nx0, ny0, simu_name[vsimu2].upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-ratio-fi-dic-alk.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 25s
2023-06-29 22:42:59.738071 Maps of the ratio of FIDIC and FIALK: plot > model: ACCESS-ESM1-5 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 >> simu: ssp126 >> simu: ssp585 > model: CESM2 >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM >> simu: picontrol >> simu: historical >> simu: ssp126 >> simu: ssp585 >> simu: ssp126 >> simu: ssp585 Figure saved: maps-ratio-fi-dic-alk.png CPU times: user 2min 25s, sys: 1.22 s, total: 2min 26s Wall time: 2min 26s
%%time
print(datetime.datetime.now())
print('Histogram')
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (18*cm2in, 20*cm2in) # width height
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# KEYWORDS DICT
#-----------------
kwhist = dict(bins=50, range=[0, 3], density=True, color='black')
kwhist2 = dict(bins=50, range=[0, 3], density=True, color='grey')
kwlines = dict(color='gray', ls='--')
#---------------------
# Plot histograms
#---------------------
for imodel, vmodel in enumerate(model_list):
for isimu, vsimu in enumerate(simu_list):
zax = ax[imodel, isimu]
if vsimu == 'historical':
zwdata = data2plot[vmodel][vsimu]['Z'].flatten()
pcm = zax.hist(zwdata, **kwhist2)
else:
zwdata = data2plot[vmodel][vsimu]['Z'].flatten()
pcm = zax.hist(zwdata, **kwhist)
zwdatahist = data2plot[vmodel]['historical']['Z'].flatten()
pcm = zax.hist(zwdatahist, alpha=.5, **kwhist2)
quantile20 = np.nanquantile(zwdata, .2)
quantile50 = np.nanquantile(zwdata, .5)
quantile80 = np.nanquantile(zwdata, .8)
zw= np.nonzero(zwdata > 1.)
pctabove50 = len(zw[0])/len(np.where(zwdata==zwdata)[0])*100
#zax.axhline(0.2, **kwlines)
#zax.axvline(quantile20, **kwlines)
zax.set_title(subnum.pop()+") %.0f%s > 1\n80%s -- 50%s -- 20%s\n%.2f -- %.2f -- %.2f"
%(pctabove50, "%", "%", "%", "%", quantile20, quantile50, quantile80), loc='left')
#
fig.tight_layout()
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .4*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .9*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'histogram-ratio-fi-dic-alk.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 7s
2022-11-21 17:31:14.871035 Histogram
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:49: RuntimeWarning: invalid value encountered in greater
Figure saved: histogram-ratio-fi-dic-alk.png CPU times: user 6.89 s, sys: 163 ms, total: 7.05 s Wall time: 6.78 s
one figure for DIC, one for ALK, one for DICsn, on for ALKsn
row = model, col = simu
%%time
print(datetime.datetime.now())
print('Maps of variance of DIC and ALK: input')
data2plot_dir = dirroot + 'dimred-220513-predictors-variances/'
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['historical', 'ssp126', 'ssp585']
var_list = ['dissic', 'dissic_sn', 'talk', 'talk_sn']
namevar_dict = {
'dissic' : 'DIC' ,
'talk' : 'ALK' ,
'dissic_sn': 'DIC-sn',
'talk_sn' : 'ALK-sn'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
zfact = 1e3 # mol.m-3 to mmol.m-3 (or umol.l)
2023-11-28 17:15:19.625043 Maps of variance of DIC and ALK: input CPU times: user 111 µs, sys: 116 µs, total: 227 µs Wall time: 205 µs
%%time
print(datetime.datetime.now())
print('Maps of variance of DIC and ALK: plot')
for vvar in var_list:
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (12*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('viridis', 256)
newcolors = zw(np.linspace(0, 1, 10+1))
cmap1 = ListedColormap(newcolors[:-1])
cmap1.set_under(color='silver')
cmap1.set_over (color=newcolors[-1])
cmap1.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'cmap':cmap1, 'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
vmin, vmax = {}, {}
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
#__________
# Load data2plot
data2plot_file = data2plot_dir + 'data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-'+vmodel+'.pckl'
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
#__________
# define vmin and vmax
zw = []
for simu in simu_list:
zw.extend(data2plot[vvar][simu]['Z'].flatten() * zfact)
#
zwvmin = 0
zwvmax = np.nanquantile(zw, 0.7)
f=np.floor(np.log10(zwvmax))
zwvmax = np.round(zwvmax*10**(-f))*10**f
vmax[vmodel] = zwvmax
vmin[vmodel] = zwvmin
#__________
# Plot
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
zax = ax[imodel, isimu]
X = data2plot[vvar][vsimu]['X']
Y = data2plot[vvar][vsimu]['Y']
Z = data2plot[vvar][vsimu]['Z'] * zfact
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, vmin = zwvmin, vmax = zwvmax, **kwmap)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
# zw1 = ax[0, 0].get_position()
# zw2 = ax[0, -1].get_position()
# nx0 = zw1.x0 + .5*zw1.width
# ny0 = zw1.y1 + 1*zw1.height
# nw = zw2.x1 - .5*zw2.width - nx0
# nh = 0.04*nw
# cax = fig.add_axes([nx0, ny0, nw, nh])
# cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
# orientation='horizontal', extend='max',
# label='Interannual standard deviation of %s [mmol$\,$m$^{-3}$]' %namevar_dict[vvar])
# cbar.set_ticks( np.linspace( kwmap['vmin'], kwmap['vmax'], np.int(mmm[0, 0].cmap.colors.shape[0]/2+1) ) )
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, -1].get_position()
nx0 = zw1.x1 + .2*zw1.width
ny0 = zw1.y0
nh = zw1.height
nw = 0.12*nh
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[imodel, -1], cax=cax, ticklocation='right',
orientation='vertical', extend='max')
cbar.set_ticks([vmin[vmodel], vmax[vmodel]])
#
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
stitle='Interannual standard deviation of %s [mmol$\,$m$^{-3}$]' %namevar_dict[vvar]
fig.suptitle(stitle, y=1.07)
#---------------------
# Save figure
#---------------------
fignam = 'maps-variance-'+namevar_dict[vvar].lower()+'.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
#
# ca. 2m30s
2022-11-21 17:32:59.948785 Maps of variance of DIC and ALK: plot > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-variance-dic.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-variance-dic-sn.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-variance-alk.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-variance-alk-sn.png CPU times: user 1min 20s, sys: 8.16 s, total: 1min 28s Wall time: 1min 7s
%%time
print(datetime.datetime.now())
print('Maps of variance of DIC and ALK: plot delta')
# var_list = ['dissic_sn', 'talk_sn']
#-----------------
# define vmin and vmax
#-----------------
vmin, vmax = {}, {}
for imodel, vmodel in enumerate(model_list):
vmin[vmodel], vmax[vmodel] = {}, {}
#__________
# Load data2plot
data2plot_file = data2plot_dir + 'data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-'+vmodel+'.pckl'
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
zwhist = []
for vvar in var_list: zwhist.extend( data2plot[vvar]['historical']['Z'].flatten() * zfact )
zwvmax = np.nanquantile(zwhist, 0.7)
f=np.floor(np.log10(zwvmax))
zwvmax = np.round(zwvmax*10**(-f))*10**f
vmax[vmodel]['hist'] = zwvmax
vmin[vmodel]['hist'] = 0
zwssp = []
for simu in ['ssp126', 'ssp585']:
for vvar in var_list:
zw = data2plot[vvar][simu]['Z'].flatten() / data2plot[vvar]['historical']['Z'].flatten() - 1
zwssp.extend( np.abs(zw) )
#
#
zwvmax = np.nanquantile(zwssp, 0.9)
f=np.floor(np.log10(zwvmax))
zwvmax = np.round(zwvmax*10**(-f))*10**f
vmax[vmodel]['ssps'] = 1+zwvmax
vmin[vmodel]['ssps'] = 1-zwvmax
#
for vvar in var_list:
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (12*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('viridis', 256)
newcolors = zw(np.linspace(0, 1, 10+1))
cmap1 = ListedColormap(newcolors[:-1])
cmap1.set_under(color='silver')
cmap1.set_over (color=newcolors[-1])
cmap1.set_bad(color='silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap2 = ListedColormap(newcolors[1:-1])
cmap2.set_under(color=newcolors[1])
cmap2.set_over (color=newcolors[-1])
cmap2.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap1 = {'cmap':cmap1, 'transform':ccrs.PlateCarree() }
kwmap2 = {'cmap':cmap2, 'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
#__________
# Load data2plot
data2plot_file = data2plot_dir + 'data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-'+vmodel+'.pckl'
with open(data2plot_file, 'rb') as f1: data2plot = pickle.load(f1)
print("File loaded: "+data2plot_file)
#__________
# Plot
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
zax = ax[imodel, isimu]
X = data2plot[vvar][vsimu]['X']
Y = data2plot[vvar][vsimu]['Y']
if vsimu == 'historical':
Z = data2plot[vvar][vsimu]['Z'] * zfact
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, vmin = vmin[vmodel]['hist'] , vmax = vmax[vmodel]['hist'], **kwmap1)
else:
Z = ( data2plot[vvar][vsimu]['Z'] / data2plot[vvar]['historical']['Z'] )
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, vmin = vmin[vmodel]['ssps'], vmax = vmax[vmodel]['ssps'], **kwmap2)
#
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
# zw1 = ax[0, 0].get_position()
# zw2 = ax[0, -1].get_position()
# nx0 = zw1.x0 + .5*zw1.width
# ny0 = zw1.y1 + 1*zw1.height
# nw = zw2.x1 - .5*zw2.width - nx0
# nh = 0.04*nw
# cax = fig.add_axes([nx0, ny0, nw, nh])
# cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
# orientation='horizontal', extend='max',
# label='Interannual standard deviation of %s [mmol$\,$m$^{-3}$]' %namevar_dict[vvar])
# cbar.set_ticks( np.linspace( kwmap['vmin'], kwmap['vmax'], np.int(mmm[0, 0].cmap.colors.shape[0]/2+1) ) )
for imodel, vmodel in enumerate(model_list):
# Hist
zw1 = ax[imodel, 0].get_position()
nx0 = zw1.x0 - .2*zw1.width
ny0 = zw1.y0
nh = zw1.height
nw = 0.12*nh
cax = fig.add_axes([nx0-nw, ny0, nw, nh])
cbar = fig.colorbar(mmm[imodel, 0], cax=cax, ticklocation='left',
orientation='vertical', extend='max')
cbar.set_ticks([vmin[vmodel]['hist'], 0, vmax[vmodel]['hist']])
# SSPs
zw1 = ax[imodel, -1].get_position()
nx0 = zw1.x1 + .2*zw1.width
ny0 = zw1.y0
nh = zw1.height
nw = 0.12*nh
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[imodel, -1], cax=cax, ticklocation='right',
orientation='vertical', extend='both')
cbar.set_ticks([vmin[vmodel]['ssps'], 1, vmax[vmodel]['ssps']])
#
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .5*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
stitle='Interannual standard deviation of %s [mmol$\,$m$^{-3}$]' %namevar_dict[vvar]
fig.suptitle(stitle, y=1.07)
#---------------------
# Save figure
#---------------------
fignam = 'maps-delta-variance-'+namevar_dict[vvar].lower()+'.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
#
# ca. 2m30s
2022-11-21 17:34:13.957546 Maps of variance of DIC and ALK: plot delta File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-delta-variance-dic.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-delta-variance-dic-sn.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-delta-variance-alk.png > model: ACCESS-ESM1-5 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: CESM2 File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: IPSL-CM6A-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: MPI-ESM1-2-LR File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 > model: NorESM2-LM File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-delta-variance-alk-sn.png CPU times: user 1min 19s, sys: 7.59 s, total: 1min 27s Wall time: 1min 7s
%%time
print(datetime.datetime.now())
print('Maps of variance of DIC and ALK: plot delta DIC vs delta ALK')
#-----------------
# Prepare data2plot
#-----------------
var_list = ['dissic_sn', 'talk_sn']
simu_list = ['ssp126', 'ssp585']
data2plot = {}
for vvar in var_list:
data2plot[vvar] = {}
for imodel, vmodel in enumerate(model_list):
data2plot[vvar][vmodel] = {}
data2plot_file = data2plot_dir + 'data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-'+vmodel+'.pckl'
with open(data2plot_file, 'rb') as f1: loadeddatas = pickle.load(f1)
print("File loaded: "+data2plot_file)
for isimu, vsimu in enumerate(simu_list):
data2plot[vvar][vmodel][vsimu] = {}
X = loadeddatas[vvar][vsimu]['X']
Y = loadeddatas[vvar][vsimu]['Y']
Z = ( loadeddatas[vvar][vsimu]['Z'] / loadeddatas[vvar]['historical']['Z'] )
data2plot[vvar][vmodel][vsimu]['X'] = X
data2plot[vvar][vmodel][vsimu]['Y'] = Y
data2plot[vvar][vmodel][vsimu]['Z'] = Z
#
#
#
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (8*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('RdBu', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap = ListedColormap(newcolors[1:-1])
cmap.set_under(color=newcolors[1])
cmap.set_over (color=newcolors[-1])
cmap.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':2, 'cmap':cmap, 'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
for isimu, vsimu in enumerate(simu_list):
zax = ax[imodel, isimu]
X = data2plot['dissic_sn'][vmodel][vsimu]['X']
Y = data2plot['dissic_sn'][vmodel][vsimu]['Y']
Z = data2plot['dissic_sn'][vmodel][vsimu]['Z']/data2plot['talk_sn'][vmodel][vsimu]['Z']
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, -1].get_position()
nx0 = zw1.x0 + .5*zw1.width
ny0 = zw1.y1 + .8*zw1.height
nw = zw2.x1 - .5*zw2.width - nx0
nh = 0.04*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='max',
label='$\Delta\,\sigma_{DIC}\,/\,\Delta\,\sigma_{ALK}$' )
cbar.set_ticks( np.linspace( kwmap['vmin'], kwmap['vmax'], np.int(mmm[0, 0].cmap.colors.shape[0]/2+1) ) )
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .3*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
stitle='$\sigma_{X}$ : inter-annual variability of X'+\
'\n$\Delta\,\sigma_{X}=\sigma_{X}^{SSP}\,/\,\sigma_{X}^{HIST}$'+\
'\nWhen < 1 (red), $\sigma_{DIC}$ decrease more (or increase less) than $\sigma_{ALK}$ in SSP'
fig.suptitle(stitle, y=1.15, va='bottom')
#---------------------
# Save figure
#---------------------
fignam = 'maps-ratio-delta-var-dic-delta-var-alk.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 13s
2022-11-21 17:35:28.245285 Maps of variance of DIC and ALK: plot delta DIC vs delta ALK File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl Figure saved: maps-ratio-delta-var-dic-delta-var-alk.png CPU times: user 14.4 s, sys: 1.63 s, total: 16 s Wall time: 11.8 s
%%time
print(datetime.datetime.now())
print('Maps of variance of DIC and ALK: plot delta DIC vs delta ALK NOT SALI NORMALIZED')
#-----------------
# Prepare data2plot
#-----------------
var_list = ['dissic', 'talk']
simu_list = ['ssp126', 'ssp585']
data2plot = {}
for vvar in var_list:
data2plot[vvar] = {}
for imodel, vmodel in enumerate(model_list):
data2plot[vvar][vmodel] = {}
data2plot_file = data2plot_dir + 'data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-'+vmodel+'.pckl'
with open(data2plot_file, 'rb') as f1: loadeddatas = pickle.load(f1)
print("File loaded: "+data2plot_file)
for isimu, vsimu in enumerate(simu_list):
data2plot[vvar][vmodel][vsimu] = {}
X = loadeddatas[vvar][vsimu]['X']
Y = loadeddatas[vvar][vsimu]['Y']
Z = ( loadeddatas[vvar][vsimu]['Z'] / loadeddatas[vvar]['historical']['Z'] )
data2plot[vvar][vmodel][vsimu]['X'] = X
data2plot[vvar][vmodel][vsimu]['Y'] = Y
data2plot[vvar][vmodel][vsimu]['Z'] = Z
#
#
#
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = len(simu_list)
fsize = (8*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('RdBu', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap = ListedColormap(newcolors[1:-1])
cmap.set_under(color=newcolors[1])
cmap.set_over (color=newcolors[-1])
cmap.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap = {'vmin':0, 'vmax':2, 'cmap':cmap, 'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
for isimu, vsimu in enumerate(simu_list):
zax = ax[imodel, isimu]
X = data2plot['dissic'][vmodel][vsimu]['X']
Y = data2plot['dissic'][vmodel][vsimu]['Y']
Z = data2plot['dissic'][vmodel][vsimu]['Z']/data2plot['talk'][vmodel][vsimu]['Z']
mmm[imodel, isimu] = zax.pcolormesh(X, Y, Z, **kwmap)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, -1].get_position()
nx0 = zw1.x0 + .5*zw1.width
ny0 = zw1.y1 + .8*zw1.height
nw = zw2.x1 - .5*zw2.width - nx0
nh = 0.04*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='max',
label='$\Delta\,\sigma_{DIC}\,/\,\Delta\,\sigma_{ALK}$' )
cbar.set_ticks( np.linspace( kwmap['vmin'], kwmap['vmax'], np.int(mmm[0, 0].cmap.colors.shape[0]/2+1) ) )
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .3*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
stitle='$\sigma_{X}$ : inter-annual variability of X'+\
'\n$\Delta\,\sigma_{X}=\sigma_{X}^{SSP}\,/\,\sigma_{X}^{HIST}$'+\
'\nWhen < 1 (red), $\sigma_{DIC}$ decrease more (or increase less) than $\sigma_{ALK}$ in SSP'
fig.suptitle(stitle, y=1.15, va='bottom')
#---------------------
# Save figure
#---------------------
fignam = 'maps-ratio-delta-var-dic-delta-var-alk-not-sali-normalized.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 13s
2023-11-28 17:15:22.816554 Maps of variance of DIC and ALK: plot delta DIC vs delta ALK NOT SALI NORMALIZED File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-ACCESS-ESM1-5.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-CESM2.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-IPSL-CM6A-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-MPI-ESM1-2-LR.pckl File loaded: /mnt/reef-ns1002k/daco/MY_JUPYTER_NOTEBOOKS/DIM-RED/dimred-220513-predictors-variances/data2plot-variances-of-dtrd-deseas-dic-alk-dicsn-alksn-NorESM2-LM.pckl
<timed exec>:119: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
Figure saved: maps-ratio-delta-var-dic-delta-var-alk-not-sali-normalized.png CPU times: user 12.9 s, sys: 919 ms, total: 13.9 s Wall time: 13.1 s
one figure for DIC, one for ALK
row = model, col = simu
%%time
print(datetime.datetime.now())
print('Maps of changes in vertical gradient of DIC and ALK: input')
data2plot_dir_dict = {'dissic': dirroot+'DIC_GRADZ/',
'talk' : dirroot+'ALK_GRADZ/'}
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
period_list = ['1985-2014', '2071-2100']
var_list = ['dissic', 'talk']
grad_depth_list = ['surf-250m', 'surf-500m']
namevar_dict = {
'dissic' : 'DIC' ,
'talk' : 'ALK'
}
2022-12-05 21:31:34.014501 Maps of changes in vertical gradient of DIC and ALK: input CPU times: user 0 ns, sys: 578 µs, total: 578 µs Wall time: 392 µs
def read_grad_dic_datas(zwfname):
zwgridshape = (360, 180)
zwdata = np.loadtxt(zwfname, dtype=str, delimiter=';')
# zwdata = np.genfromtxt(zwfname, dtype=str)
Z = np.array(np.where(zwdata.flatten()=='NA', np.nan, zwdata.flatten()).reshape(zwgridshape), dtype=float).T
# X and Y are adapted for plt.pcolormesh with shading='flat'
zwX = np.arange(-179.5, 179.5+1/2., 1) + 179.5
zwY = np.arange(-89.5, 89.5+1/2., 1)
X, Y = np.meshgrid(zwX, zwY)
return Z, Y, X
%%time
print(datetime.datetime.now())
print('Maps of changes in vertical gradient of DIC and ALK: plot')
for vvar in var_list:
for vgrad_depth in grad_depth_list:
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = 3
fsize = (12*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP AND KEYWORDS DICT
#-----------------
if vvar=='dissic':
#_________
# cmap
nbcolors = 9
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors+1+4:]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors]
cmap1 = ListedColormap(aa[:-1])
cmap1.set_under(color = bb)
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+2))[:nbcolors+1]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+2))[nbcolors+1+2]
cmap2 = ListedColormap(aa[1:])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = bb)
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1 = {'vmin':0, 'vmax':.9, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2 = {'vmin':-.6, 'vmax':0, 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1 = np.linspace( kwmap1['vmin'], kwmap1['vmax'], np.int(cmap1.colors.shape[0]/3+1) )
cbticks2 = np.linspace( kwmap2['vmin'], kwmap2['vmax'], np.int(cmap2.colors.shape[0]/3+1) )
elif vvar=='talk':
#_________
# cmap
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 10))
cmap1 = ListedColormap(aa[1:-1])
cmap1.set_under(color = aa[0])
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 10))
cmap2 = ListedColormap(aa[1:-1])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = aa[-1])
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1 = {'vmin':-.4, 'vmax':.4, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2 = {'vmin':-.2, 'vmax':.2, 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1 = np.linspace( kwmap1['vmin'], kwmap1['vmax'], np.int(cmap1.colors.shape[0]/2+1) )
cbticks2 = np.linspace( kwmap2['vmin'], kwmap2['vmax'], np.int(cmap2.colors.shape[0]/2+1) )
else: exit('Wrong vvar')
kwcont = dict(levels=[0], colors='k', linewidths=1., linestyles='--', transform=ccrs.PlateCarree())
#---------------------
# MAPS
#---------------------
pcm = np.zeros_like(ax)
for imodel, vmodel in enumerate(model_list):
zwsaved = {}
for iperiod, vperiod in enumerate(period_list):
#____________
# Load datas
fname = data2plot_dir_dict[vvar]+vvar+'_'+vmodel+'_'+vperiod+'_Grad'+vgrad_depth+'.txt'
zw = read_grad_dic_datas(fname)
if vgrad_depth == 'surf-500m' : ddd = 500.
elif vgrad_depth == 'surf-250m' : ddd = 250.
else: exit('grad_depth not correct')
X = zw[2]
Y = zw[1]
Z = -zw[0]/ddd
zwsaved[vperiod] = Z
#____________
# Plot
zax = ax[imodel, iperiod]
pcm[imodel, iperiod] = zax.pcolormesh(X, Y, Z, **kwmap1)
zax.contour(X, Y, Z, **kwcont)
zax.coastlines()
#
#____________
# Plot delta
zax = ax[imodel, iperiod+1]
Z = zwsaved[period_list[1]] - zwsaved[period_list[0]]
pcm[imodel, iperiod+1] = zax.pcolormesh(X, Y, Z, **kwmap2)
zax.contour(X, Y, Z, **kwcont)
zax.coastlines()
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
zw1 = ax[0, 0].get_position()
zw2 = ax[0, 1].get_position()
nx0 = zw1.x0 + .5*zw1.width
ny0 = zw1.y1 + 1*zw1.height
nw = zw2.x1 - .5*zw2.width - nx0
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
zwlab = '%s vertical gradient (%s)\n in [mmol$\,$m$^{-3}\,$m$^{-1}$]' %(namevar_dict[vvar], vgrad_depth)
cbar = fig.colorbar(pcm[0, 0], cax=cax, ticklocation='top',
orientation='horizontal', extend='both',
label = zwlab)
cbar.set_ticks( cbticks1 )
zw1 = ax[0, -1].get_position()
nx0 = zw1.x0
ny0 = zw1.y1 + 1*zw1.height
nw = zw1.width
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
zwlab = '[mmol$\,$m$^{-3}\,$m$^{-1}$]'
cbar = fig.colorbar(pcm[0, -1], cax=cax, ticklocation='top',
orientation='horizontal', extend='both',
label = zwlab)
cbar.set_ticks( cbticks2 )
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: period names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for iperiod, vperiod in enumerate(period_list):
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
#
iperiod, vperiod = 2, 'Change'
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
#---------------------
# Save figure
#---------------------
fignam = 'maps-gradz-'+namevar_dict[vvar].lower()+'-'+vgrad_depth.lower()+'.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
#
#
# ca. 1m
2022-11-21 17:36:32.948900 Maps of changes in vertical gradient of DIC and ALK: plot Figure saved: maps-gradz-dic-surf-250m.png
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1366: UserWarning: No contour levels were found within the data range. result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
Figure saved: maps-gradz-dic-surf-500m.png Figure saved: maps-gradz-alk-surf-250m.png Figure saved: maps-gradz-alk-surf-500m.png CPU times: user 1min 21s, sys: 9.94 s, total: 1min 31s Wall time: 1min 6s
%%time
print(datetime.datetime.now())
print('Maps of changes in vertical gradient of DIC and ALK: plot summary')
vgrad_depth = 'surf-500m'
gradname = '500m-surf'
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = 4
fsize = (16*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP AND KEYWORDS DICT
#-----------------
kwmap1, kwmap2 = {}, {}
cbticks1, cbticks2 = {}, {}
for ivar, vvar in enumerate(var_list):
if vvar=='dissic':
#_________
# cmap
nbcolors = 9
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors+1+4:]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors]
cmap1 = ListedColormap(aa[:-1])
cmap1.set_under(color = bb)
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+2))[:nbcolors+1]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+2))[nbcolors+1+2]
cmap2 = ListedColormap(aa[1:])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = bb)
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1[vvar] = {'vmin':0, 'vmax':.9, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2[vvar] = {'vmin':-.6, 'vmax':0, 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1[vvar] = np.linspace( kwmap1[vvar]['vmin'], kwmap1[vvar]['vmax'], np.int(cmap1.colors.shape[0]/3+1) )
cbticks2[vvar] = np.linspace( kwmap2[vvar]['vmin'], kwmap2[vvar]['vmax'], np.int(cmap2.colors.shape[0]/3+1) )
elif vvar=='talk':
#_________
# cmap
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 10))
cmap1 = ListedColormap(aa[1:-1])
cmap1.set_under(color = aa[0])
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 10))
cmap2 = ListedColormap(aa[1:-1])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = aa[-1])
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1[vvar] = {'vmin':-.4, 'vmax':.4, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2[vvar] = {'vmin':-.2, 'vmax':.2, 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1[vvar] = np.linspace( kwmap1[vvar]['vmin'], kwmap1[vvar]['vmax'], np.int(cmap1.colors.shape[0]/2+1) )
cbticks2[vvar] = np.linspace( kwmap2[vvar]['vmin'], kwmap2[vvar]['vmax'], np.int(cmap2.colors.shape[0]/2+1) )
else: exit('Wrong vvar')
#
kwcont = dict(levels=[0], colors='k', linewidths=1., linestyles='--', transform=ccrs.PlateCarree())
#---------------------
# MAPS
#---------------------
pcm = np.zeros_like(ax)
for ivar, vvar in enumerate(var_list):
for imodel, vmodel in enumerate(model_list):
if vgrad_depth == 'surf-500m' : ddd = 500.
elif vgrad_depth == 'surf-250m' : ddd = 250.
else: exit('grad_depth not correct')
#____________
# Load datas
vperiod = '1985-2014'
fname = data2plot_dir_dict[vvar]+vvar+'_'+vmodel+'_'+vperiod+'_Grad'+vgrad_depth+'.txt'
zw = read_grad_dic_datas(fname)
X = zw[2]
Y = zw[1]
Zhist = -zw[0]/ddd
vperiod = '2071-2100'
fname = data2plot_dir_dict[vvar]+vvar+'_'+vmodel+'_'+vperiod+'_Grad'+vgrad_depth+'.txt'
zw = read_grad_dic_datas(fname)
Zchange = -zw[0]/ddd - Zhist
#____________
# Plot
iperiod = 0 + 2*ivar
zax = ax[imodel, iperiod]
pcm[imodel, iperiod] = zax.pcolormesh(X, Y, Zhist, **kwmap1[vvar])
zax.contour(X, Y, Zhist, **kwcont)
zax.coastlines()
iperiod = 1 + 2*ivar
zax = ax[imodel, iperiod]
pcm[imodel, iperiod] = zax.pcolormesh(X, Y, Zchange, **kwmap2[vvar])
zax.contour(X, Y, Zchange, **kwcont)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
for ivar, vvar in enumerate(var_list):
iperiod = 2*ivar
zw1 = ax[-1, iperiod].get_position()
nx0 = zw1.x0 + .1*zw1.width
ny0 = zw1.y0 - .3*zw1.height
nw = .8*zw1.width
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm[0, iperiod], cax=cax, ticklocation='bottom',
orientation='horizontal', extend='both')
cbar.set_ticks( cbticks1[vvar] )
iperiod = 2*ivar+1
zw1 = ax[-1, iperiod].get_position()
nx0 = zw1.x0 + .1*zw1.width
ny0 = zw1.y0 - .3*zw1.height
nw = .8*zw1.width
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm[0, iperiod], cax=cax, ticklocation='bottom',
orientation='horizontal', extend='both')
cbar.set_ticks( cbticks2[vvar] )
iperiod = 2*ivar
zw1 = ax[-1, iperiod].get_position()
iperiod = 2*ivar+1
zw2 = ax[-1, iperiod].get_position()
nx0 = 0.5*(zw1.x1 + zw2.x0)
ny0 = zw1.y0 - 1.3*zw1.height
kwtxt = dict(fontsize='large', rotation='horizontal', ha='center', va='bottom')
zwlab = '%s vertical gradient (%s)\n in [mmol$\,$m$^{-3}\,$m$^{-1}$]' %(namevar_dict[vvar], gradname)
fig.text(nx0, ny0, zwlab, **kwtxt)
#
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: period names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for ivar, vvar in enumerate(var_list):
iperiod, vperiod = 2*ivar, period_list[0]
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
iperiod, vperiod = 2*ivar+1, 'Change'
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-gradz-alk-dic-summary.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 1m
2022-12-05 21:32:44.604610 Maps of changes in vertical gradient of DIC and ALK: plot summary Figure saved: maps-gradz-alk-dic-summary.png CPU times: user 24.7 s, sys: 2.6 s, total: 27.3 s Wall time: 20.1 s
%%time
print(datetime.datetime.now())
print('Maps of changes in vertical gradient of DIC and ALK: plot summary 2')
vgrad_depth = 'surf-500m'
gradname = '500m-surf'
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(model_list)
ncol = 4
fsize = (16*cm2in, 14*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP AND KEYWORDS DICT
#-----------------
kwmap1, kwmap2 = {}, {}
cbticks1, cbticks2 = {}, {}
for ivar, vvar in enumerate(var_list):
if vvar=='dissic':
#_________
# cmap
nbcolors = 9
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors+1+4:]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+4))[nbcolors]
cmap1 = ListedColormap(aa[:-1])
cmap1.set_under(color = bb)
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 2*nbcolors+2+2))[:nbcolors+1]
bb = zw(np.linspace(0, 1, 2*nbcolors+2+2))[nbcolors+1+2]
cmap2 = ListedColormap(aa[1:])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = bb)
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1[vvar] = {'vmin':0, 'vmax':.9, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2[vvar] = {'vmin':-1., 'vmax':0., 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1[vvar] = np.linspace( kwmap1[vvar]['vmin'], kwmap1[vvar]['vmax'], np.int(cmap1.colors.shape[0]/3+1) )
cbticks2[vvar] = np.linspace( kwmap2[vvar]['vmin'], kwmap2[vvar]['vmax'], np.int(cmap2.colors.shape[0]/3+1) )
elif vvar=='talk':
#_________
# cmap
zw = cm.get_cmap('PiYG', 256)
aa = zw(np.linspace(0, 1, 10))
cmap1 = ListedColormap(aa[1:-1])
cmap1.set_under(color = aa[0])
cmap1.set_over(color = aa[-1])
cmap1.set_bad(color = 'silver', alpha=0)
zw = cm.get_cmap('RdBu', 256)
aa = zw(np.linspace(0, 1, 10))
cmap2 = ListedColormap(aa[1:-1])
cmap2.set_under(color = aa[0])
cmap2.set_over(color = aa[-1])
cmap2.set_bad(color = 'silver', alpha=0)
#_________
# kwmap
kwmap1[vvar] = {'vmin':-.4, 'vmax':.4, 'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2[vvar] = {'vmin':-1, 'vmax':1, 'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
cbticks1[vvar] = np.linspace( kwmap1[vvar]['vmin'], kwmap1[vvar]['vmax'], np.int(cmap1.colors.shape[0]/2+1) )
cbticks2[vvar] = np.linspace( kwmap2[vvar]['vmin'], kwmap2[vvar]['vmax'], np.int(cmap2.colors.shape[0]/2+1) )
else: exit('Wrong vvar')
#
kwcont = dict(levels=[0], colors='k', linewidths=1., linestyles='--', transform=ccrs.PlateCarree())
#---------------------
# MAPS
#---------------------
pcm = np.zeros_like(ax)
for ivar, vvar in enumerate(var_list):
for imodel, vmodel in enumerate(model_list):
if vgrad_depth == 'surf-500m' : ddd = 500.
elif vgrad_depth == 'surf-250m' : ddd = 250.
else: exit('grad_depth not correct')
#____________
# Load datas
vperiod = '1985-2014'
fname = data2plot_dir_dict[vvar]+vvar+'_'+vmodel+'_'+vperiod+'_Grad'+vgrad_depth+'.txt'
zw = read_grad_dic_datas(fname)
X = zw[2]
Y = zw[1]
Zhist = -zw[0]/ddd
vperiod = '2071-2100'
fname = data2plot_dir_dict[vvar]+vvar+'_'+vmodel+'_'+vperiod+'_Grad'+vgrad_depth+'.txt'
zw = read_grad_dic_datas(fname)
Zchange = (-zw[0]/ddd - Zhist)/Zhist
#____________
# Plot
iperiod = 0 + 2*ivar
zax = ax[imodel, iperiod]
pcm[imodel, iperiod] = zax.pcolormesh(X, Y, Zhist, **kwmap1[vvar])
zax.contour(X, Y, Zhist, **kwcont)
zax.coastlines()
iperiod = 1 + 2*ivar
zax = ax[imodel, iperiod]
pcm[imodel, iperiod] = zax.pcolormesh(X, Y, Zchange, **kwmap2[vvar])
zax.contour(X, Y, Zchange, **kwcont)
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for icol, zax in enumerate(axrow):
zwtitle = subnum.pop() + ')'
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
for ivar, vvar in enumerate(var_list):
iperiod = 2*ivar
zw1 = ax[0, iperiod].get_position()
nx0 = zw1.x0 + .1*zw1.width
ny0 = zw1.y1 + 1*zw1.height
nw = .8*zw1.width
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm[0, iperiod], cax=cax, ticklocation='top',
orientation='horizontal', extend='both')
cbar.set_ticks( cbticks1[vvar] )
iperiod = 2*ivar+1
zw1 = ax[0, iperiod].get_position()
nx0 = zw1.x0 + .1*zw1.width
ny0 = zw1.y1 + 1*zw1.height
nw = .8*zw1.width
nh = 0.06*nw
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(pcm[0, iperiod], cax=cax, ticklocation='top',
orientation='horizontal', extend='both')
cbar.set_ticks( cbticks2[vvar] )
iperiod = 2*ivar
zw1 = ax[0, iperiod].get_position()
iperiod = 2*ivar+1
zw2 = ax[0, iperiod].get_position()
nx0 = 0.5*(zw1.x1 + zw2.x0)
ny0 = zw1.y1 + 1.5*zw1.height
kwtxt = dict(fontsize='large', rotation='horizontal', ha='center', va='bottom')
zwlab = '%s vertical gradient (%s)\n in [mmol$\,$m$^{-3}\,$m$^{-1}$]' %(namevar_dict[vvar], gradname)
fig.text(nx0, ny0, zwlab, **kwtxt)
#
#---------------------
# Add text on the left: model names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for imodel, vmodel in enumerate(model_list):
zw1 = ax[imodel, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, vmodel.upper(), **kwtxt)
#
#---------------------
# Add text on the top: period names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for ivar, vvar in enumerate(var_list):
iperiod, vperiod = 2*ivar, period_list[0]
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
iperiod, vperiod = 2*ivar+1, 'Change'
zw1 = ax[0, iperiod].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, vperiod.upper(), **kwtxt)
#
#---------------------
# Save figure
#---------------------
fignam = 'maps-gradz-alk-dic-summary-2.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
# ca. 1M
2022-09-21 15:42:18.916015 Maps of changes in vertical gradient of DIC and ALK: plot summary 2
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:125: RuntimeWarning: divide by zero encountered in true_divide
Figure saved: maps-gradz-alk-dic-summary-2.png CPU times: user 30 s, sys: 4.09 s, total: 34.1 s Wall time: 24.9 s
%%time
print(datetime.datetime.now())
print('Maps of predictors and cflux: inputs')
model_list = ['ACCESS-ESM1-5', 'CESM2', 'IPSL-CM6A-LR', 'MPI-ESM1-2-LR', 'NorESM2-LM']
simu_list = ['historical', 'ssp126', 'ssp585']
predictor_list = ['fgco2', 'dissic', 'talk', 'tos', 'sfcWind', 'siconc', 'sos', 'po4']
feature_name = {
'siconc' :'siconc',
'talk' :'ALK',
'dissic' :'DIC',
'tos' :'SST',
'sfcWind':'sfcWind',
'sos' :'SSS',
'po4' :'PO4',
'atmco2' :'atmCO2',
'fgco2' :'CO$_2$ flux'
}
simu_name = {
'picontrol': 'preindustrial',
'historical': 'historical',
'ssp126': 'SSP126',
'ssp585': 'SSP585'
}
2022-11-21 17:45:56.544424 Maps of predictors and cflux: inputs CPU times: user 0 ns, sys: 587 µs, total: 587 µs Wall time: 398 µs
%%time
print(datetime.datetime.now())
print('Maps of predictors and cflux: prepare data2plot')
data2plot = {}
for model in model_list:
print('> model: '+model)
data2plot[model]={}
for simu in simu_list:
print('>> simu: '+simu)
data2plot[model][simu] = {}
for predictor in predictor_list:
data2plot[model][simu][predictor] = {}
dirname = dircmip6 + model + '/' + simu + '/'
fname = predictor + '_' + model + '_' + simu + '*.nc'
dataset = xr.open_mfdataset(dirname + fname, use_cftime=True)
print('File loaded: %s' %fname)
zwpred = dataset[predictor].sel(time=slice('1850', '2100')).mean(dim='time').squeeze()
zwX = zwpred['lon'].values
zwY = zwpred['lat'].values
zwZ = zwpred.values
X, Y = np.meshgrid(zwX, zwY)
data2plot[model][simu][predictor]['X'] = X
data2plot[model][simu][predictor]['Y'] = Y
data2plot[model][simu][predictor]['Z'] = zwZ
#
#
#
# ca. 2m30s
2022-11-21 17:43:08.321342 Maps of predictors and cflux: prepare data2plot > model: ACCESS-ESM1-5 >> simu: historical File loaded: fgco2_ACCESS-ESM1-5_historical*.nc File loaded: dissic_ACCESS-ESM1-5_historical*.nc File loaded: talk_ACCESS-ESM1-5_historical*.nc File loaded: tos_ACCESS-ESM1-5_historical*.nc File loaded: sfcWind_ACCESS-ESM1-5_historical*.nc File loaded: siconc_ACCESS-ESM1-5_historical*.nc File loaded: sos_ACCESS-ESM1-5_historical*.nc File loaded: po4_ACCESS-ESM1-5_historical*.nc >> simu: ssp126 File loaded: fgco2_ACCESS-ESM1-5_ssp126*.nc File loaded: dissic_ACCESS-ESM1-5_ssp126*.nc File loaded: talk_ACCESS-ESM1-5_ssp126*.nc File loaded: tos_ACCESS-ESM1-5_ssp126*.nc File loaded: sfcWind_ACCESS-ESM1-5_ssp126*.nc File loaded: siconc_ACCESS-ESM1-5_ssp126*.nc File loaded: sos_ACCESS-ESM1-5_ssp126*.nc File loaded: po4_ACCESS-ESM1-5_ssp126*.nc >> simu: ssp585 File loaded: fgco2_ACCESS-ESM1-5_ssp585*.nc File loaded: dissic_ACCESS-ESM1-5_ssp585*.nc File loaded: talk_ACCESS-ESM1-5_ssp585*.nc File loaded: tos_ACCESS-ESM1-5_ssp585*.nc File loaded: sfcWind_ACCESS-ESM1-5_ssp585*.nc File loaded: siconc_ACCESS-ESM1-5_ssp585*.nc File loaded: sos_ACCESS-ESM1-5_ssp585*.nc File loaded: po4_ACCESS-ESM1-5_ssp585*.nc > model: CESM2 >> simu: historical File loaded: fgco2_CESM2_historical*.nc File loaded: dissic_CESM2_historical*.nc File loaded: talk_CESM2_historical*.nc File loaded: tos_CESM2_historical*.nc File loaded: sfcWind_CESM2_historical*.nc File loaded: siconc_CESM2_historical*.nc File loaded: sos_CESM2_historical*.nc File loaded: po4_CESM2_historical*.nc >> simu: ssp126 File loaded: fgco2_CESM2_ssp126*.nc File loaded: dissic_CESM2_ssp126*.nc File loaded: talk_CESM2_ssp126*.nc File loaded: tos_CESM2_ssp126*.nc File loaded: sfcWind_CESM2_ssp126*.nc File loaded: siconc_CESM2_ssp126*.nc File loaded: sos_CESM2_ssp126*.nc File loaded: po4_CESM2_ssp126*.nc >> simu: ssp585 File loaded: fgco2_CESM2_ssp585*.nc File loaded: dissic_CESM2_ssp585*.nc File loaded: talk_CESM2_ssp585*.nc File loaded: tos_CESM2_ssp585*.nc File loaded: sfcWind_CESM2_ssp585*.nc File loaded: siconc_CESM2_ssp585*.nc File loaded: sos_CESM2_ssp585*.nc File loaded: po4_CESM2_ssp585*.nc > model: IPSL-CM6A-LR >> simu: historical File loaded: fgco2_IPSL-CM6A-LR_historical*.nc File loaded: dissic_IPSL-CM6A-LR_historical*.nc File loaded: talk_IPSL-CM6A-LR_historical*.nc File loaded: tos_IPSL-CM6A-LR_historical*.nc File loaded: sfcWind_IPSL-CM6A-LR_historical*.nc File loaded: siconc_IPSL-CM6A-LR_historical*.nc File loaded: sos_IPSL-CM6A-LR_historical*.nc File loaded: po4_IPSL-CM6A-LR_historical*.nc >> simu: ssp126 File loaded: fgco2_IPSL-CM6A-LR_ssp126*.nc File loaded: dissic_IPSL-CM6A-LR_ssp126*.nc File loaded: talk_IPSL-CM6A-LR_ssp126*.nc File loaded: tos_IPSL-CM6A-LR_ssp126*.nc File loaded: sfcWind_IPSL-CM6A-LR_ssp126*.nc File loaded: siconc_IPSL-CM6A-LR_ssp126*.nc File loaded: sos_IPSL-CM6A-LR_ssp126*.nc File loaded: po4_IPSL-CM6A-LR_ssp126*.nc >> simu: ssp585 File loaded: fgco2_IPSL-CM6A-LR_ssp585*.nc File loaded: dissic_IPSL-CM6A-LR_ssp585*.nc File loaded: talk_IPSL-CM6A-LR_ssp585*.nc File loaded: tos_IPSL-CM6A-LR_ssp585*.nc File loaded: sfcWind_IPSL-CM6A-LR_ssp585*.nc File loaded: siconc_IPSL-CM6A-LR_ssp585*.nc File loaded: sos_IPSL-CM6A-LR_ssp585*.nc File loaded: po4_IPSL-CM6A-LR_ssp585*.nc > model: MPI-ESM1-2-LR >> simu: historical File loaded: fgco2_MPI-ESM1-2-LR_historical*.nc File loaded: dissic_MPI-ESM1-2-LR_historical*.nc File loaded: talk_MPI-ESM1-2-LR_historical*.nc File loaded: tos_MPI-ESM1-2-LR_historical*.nc File loaded: sfcWind_MPI-ESM1-2-LR_historical*.nc File loaded: siconc_MPI-ESM1-2-LR_historical*.nc File loaded: sos_MPI-ESM1-2-LR_historical*.nc File loaded: po4_MPI-ESM1-2-LR_historical*.nc >> simu: ssp126 File loaded: fgco2_MPI-ESM1-2-LR_ssp126*.nc File loaded: dissic_MPI-ESM1-2-LR_ssp126*.nc File loaded: talk_MPI-ESM1-2-LR_ssp126*.nc File loaded: tos_MPI-ESM1-2-LR_ssp126*.nc File loaded: sfcWind_MPI-ESM1-2-LR_ssp126*.nc File loaded: siconc_MPI-ESM1-2-LR_ssp126*.nc File loaded: sos_MPI-ESM1-2-LR_ssp126*.nc File loaded: po4_MPI-ESM1-2-LR_ssp126*.nc >> simu: ssp585 File loaded: fgco2_MPI-ESM1-2-LR_ssp585*.nc File loaded: dissic_MPI-ESM1-2-LR_ssp585*.nc File loaded: talk_MPI-ESM1-2-LR_ssp585*.nc File loaded: tos_MPI-ESM1-2-LR_ssp585*.nc File loaded: sfcWind_MPI-ESM1-2-LR_ssp585*.nc File loaded: siconc_MPI-ESM1-2-LR_ssp585*.nc File loaded: sos_MPI-ESM1-2-LR_ssp585*.nc File loaded: po4_MPI-ESM1-2-LR_ssp585*.nc > model: NorESM2-LM >> simu: historical File loaded: fgco2_NorESM2-LM_historical*.nc File loaded: dissic_NorESM2-LM_historical*.nc File loaded: talk_NorESM2-LM_historical*.nc File loaded: tos_NorESM2-LM_historical*.nc File loaded: sfcWind_NorESM2-LM_historical*.nc File loaded: siconc_NorESM2-LM_historical*.nc File loaded: sos_NorESM2-LM_historical*.nc File loaded: po4_NorESM2-LM_historical*.nc >> simu: ssp126 File loaded: fgco2_NorESM2-LM_ssp126*.nc File loaded: dissic_NorESM2-LM_ssp126*.nc File loaded: talk_NorESM2-LM_ssp126*.nc File loaded: tos_NorESM2-LM_ssp126*.nc File loaded: sfcWind_NorESM2-LM_ssp126*.nc File loaded: siconc_NorESM2-LM_ssp126*.nc File loaded: sos_NorESM2-LM_ssp126*.nc File loaded: po4_NorESM2-LM_ssp126*.nc >> simu: ssp585 File loaded: fgco2_NorESM2-LM_ssp585*.nc File loaded: dissic_NorESM2-LM_ssp585*.nc File loaded: talk_NorESM2-LM_ssp585*.nc File loaded: tos_NorESM2-LM_ssp585*.nc File loaded: sfcWind_NorESM2-LM_ssp585*.nc File loaded: siconc_NorESM2-LM_ssp585*.nc File loaded: sos_NorESM2-LM_ssp585*.nc File loaded: po4_NorESM2-LM_ssp585*.nc CPU times: user 1min 20s, sys: 40.7 s, total: 2min 1s Wall time: 2min 48s
%%time
print(datetime.datetime.now())
print('Maps of predictors and cflux: plot')
# define vmin and vmax0
vmin, vmax = {}, {}
for pred in predictor_list:
zw = []
for model in model_list:
for simu in simu_list:
zw.extend(data2plot[model][simu][pred]['Z'].flatten())
#
#
zwvmin = np.nanquantile(zw, 0.1)
zwvmax = np.nanquantile(zw, 0.9)
if (zwvmin>=0) & (zwvmax>=0):
ddd= zwvmax - zwvmin
f=np.floor(np.log10(ddd))
zwvmin = np.round(zwvmin*10**(-f))*10**f
zwvmax = np.round(zwvmax*10**(-f))*10**f
ddd= zwvmax - zwvmin
if zwvmin < 0.1*ddd : zwvmin=0
elif (zwvmin<0) & (zwvmax>=0):
zw = np.nanmax([np.abs(zwvmin), np.abs(zwvmax)])
f=np.floor(np.log10(zw))
zw = np.round(zw*10**(-f))*10**f
zwvmin = -zw
zwvmax=zw
else: exit('check vmin and vmax')
vmin[pred] = zwvmin
vmax[pred] = zwvmax
#
for imodel, vmodel in enumerate(model_list):
print('> model: '+vmodel)
#-----------------
# FIGURE PARAM
#-----------------
nrow = len(predictor_list)
ncol = len(simu_list)
fsize = (14*cm2in, 23*cm2in) # width height
ccrsproj = ccrs.Robinson(central_longitude=-155)
fig, ax = plt.subplots(nrow, ncol, figsize=fsize,
subplot_kw=dict(projection=ccrsproj),
squeeze = False)
subnum = list('abcdefghijklmnopqrstuvwxyz')
aaa=[]
for sss in subnum: aaa.append('a'+sss)
subnum=subnum+aaa
subnum.reverse()
#-----------------
# CREATE CUSTOM CMAP
#-----------------
zw = cm.get_cmap('viridis', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap1 = ListedColormap(newcolors[1:-1])
cmap1.set_under(color=newcolors[0])
cmap1.set_over (color=newcolors[-1])
cmap1.set_bad(color='silver', alpha=0)
zw = cm.get_cmap('RdBu_r', 256)
newcolors = zw(np.linspace(0, 1, 10+2))
cmap2 = ListedColormap(newcolors[1:-1])
cmap2.set_under(color=newcolors[0])
cmap2.set_over (color=newcolors[-1])
cmap2.set_bad(color='silver', alpha=0)
#-----------------
# KEYWORDS DICT
#-----------------
kwmap1 = {'cmap':cmap1, \
'transform':ccrs.PlateCarree() }
kwmap2 = {'cmap':cmap2, \
'transform':ccrs.PlateCarree() }
#-----------------
# Maps
#-----------------
mmm = np.zeros_like(ax)
for isimu, vsimu in enumerate(simu_list):
print('>> simu: ' + vsimu)
for ipred, vpred in enumerate(predictor_list):
zax = ax[ipred, isimu]
X = data2plot[vmodel][vsimu][vpred]['X']
Y = data2plot[vmodel][vsimu][vpred]['Y']
Z = data2plot[vmodel][vsimu][vpred]['Z']
zwvmin, zwvmax = vmin[vpred], vmax[vpred]
if (zwvmin>=0) & (zwvmax>=0):
mmm[ipred, isimu] = zax.pcolormesh(X, Y, Z, vmin=zwvmin, vmax=zwvmax, **kwmap1)
elif (zwvmin<0) & (zwvmax>=0):
mmm[ipred, isimu] = zax.pcolormesh(X, Y, Z, vmin=zwvmin, vmax=zwvmax, **kwmap2)
else: exit('check vmin and vmax')
zax.coastlines()
#
#
fig.tight_layout()
#---------------------
# Add numbering to titles of subplots
#---------------------
for irow, axrow in enumerate(ax):
for zax in axrow:
zwtitle = subnum.pop() + ') '
zax.set_title(zwtitle, loc='left')
#
#
#---------------------
# Colorbar
#---------------------
for ipred, vpred in enumerate(predictor_list):
zw1 = ax[ipred, -1].get_position()
nx0 = zw1.x1 + .2*zw1.width
ny0 = zw1.y0
nh = zw1.height
nw = 0.12*nh
cax = fig.add_axes([nx0, ny0, nw, nh])
cbar = fig.colorbar(mmm[ipred, -1], cax=cax, ticklocation='right',
orientation='vertical', extend='both')
zwvmin, zwvmax = vmin[vpred], vmax[vpred]
if (zwvmin>=0) & (zwvmax>=0):
cbar.set_ticks([zwvmin, zwvmax])
elif (zwvmin<0) & (zwvmax>=0):
cbar.set_ticks([zwvmin, 0, zwvmax])
else: exit('check vmin and vmax')
#---------------------
# Add text on the left: feature names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='right', va='center')
for ipred, vpred in enumerate(predictor_list):
zw1 = ax[ipred, 0].get_position()
ny0 = 0.5 *(zw1.y0+zw1.y1)
nx0 = zw1.x0 - .2*zw1.width
fig.text(nx0, ny0, feature_name[vpred].upper(), **kwtxt)
#
#---------------------
# Add text on the top: simu names
#---------------------
kwtxt = dict(fontsize='large', rotation='horizontal', weight='bold', ha='center', va='bottom')
for isimu, vsimu in enumerate(simu_list):
zw1 = ax[0, isimu].get_position()
nx0 = 0.5 *(zw1.x0+zw1.x1)
ny0 = zw1.y1 + .4*zw1.height
fig.text(nx0, ny0, simu_name[vsimu].upper(), **kwtxt)
#
fig.suptitle(vmodel, y=1.07)
#---------------------
# Save figure
#---------------------
fignam = 'maps-predictors-'+vmodel+'.png'
fig.savefig(dirout+fignam, bbox_inches='tight')
print('Figure saved: '+fignam)
#
# ca. 2m30s
2022-11-21 17:45:56.572609 Maps of predictors and cflux: plot > model: ACCESS-ESM1-5 >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-predictors-ACCESS-ESM1-5.png > model: CESM2 >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-predictors-CESM2.png > model: IPSL-CM6A-LR >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-predictors-IPSL-CM6A-LR.png > model: MPI-ESM1-2-LR >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-predictors-MPI-ESM1-2-LR.png > model: NorESM2-LM >> simu: historical >> simu: ssp126 >> simu: ssp585 Figure saved: maps-predictors-NorESM2-LM.png CPU times: user 3min, sys: 16.5 s, total: 3min 16s Wall time: 2min 33s